/*
 * Alias .obj loader for processing
 * programmed by Tatsuya SAITO / UCLA Design | Media Arts 
 * Created on 2005/04/17
 *
 * 
 *
 */

package npr;

import processing.core.*;
import java.awt.event.KeyEvent;
import java.awt.event.MouseEvent;
import java.io.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.net.MalformedURLException;
import java.net.URL;
import java.util.Hashtable;
import java.util.Vector;

import javax.media.opengl.*;

import npr.NPR;
import java.lang.Math;
import java.util.Collections;

/**
 * 
 * @author tatsuyas
 * @author Katrin Lang
 * 
 * TODO:
 * 
 */

public class OBJModel {

    // global variables
    Vector<PVector> vertices; // vertices

    Vector<PVector> texturev; // texture coordinates

    Vector<PVector> normv;

    Hashtable materials;

    Hashtable groups;

    Vector<ModelSegment> modelSegments;

    String objName = "default";

    String defaultMaterialName = "default";

    Group defaultGroup = new Group("default");

    ModelSegment defaultModelSegment = new ModelSegment();

    // processing variables
    static PApplet parent;

    PImage texture; // texture image applied from the code.

    // runtime rendering variables
    int mode = PApplet.POLYGON; // render mode (ex. POLYGON, POINTS ..)

    boolean flagTexture = true;

    Debug debug;
    
    private NPR renderer;

    PVector minBounds, maxBounds;

    int displayList= -1, boxyRidges= -1, boxyValleys= -1, boxyRidgesHidden= -1, boxyValleysHidden= -1, smoothRidges= -1, smoothValleys= -1, smoothRidgesHidden= -1, smoothValleysHidden= -1, apparentRidges= -1;

    GL gl;

    float boxyRidgesThreshold= -1.0f;
    float boxyRidgesRed= -1.0f, boxyRidgesGreen= -1.0f, boxyRidgesBlue= -1.0f, boxyRidgesAlpha= -1.0f;

    float boxyValleysThreshold= -1.0f;
    float boxyValleysRed= -1.0f, boxyValleysGreen= -1.0f, boxyValleysBlue= -1.0f, boxyValleysAlpha= -1.0f;

    float boxyRidgesHiddenThreshold= -1.0f;
    float boxyRidgesHiddenRed= -1.0f, boxyRidgesHiddenGreen= -1.0f, boxyRidgesHiddenBlue= -1.0f, boxyRidgesHiddenAlpha= -1.0f;

    float boxyValleysHiddenThreshold= -1.0f;
    float boxyValleysHiddenRed= -1.0f, boxyValleysHiddenGreen= -1.0f, boxyValleysHiddenBlue= -1.0f, boxyValleysHiddenAlpha= -1.0f;


    boolean smoothRidgesDrawFaded;
    float smoothRidgesThreshold= -1.0f, smoothRidgesRed= -1.0f, smoothRidgesGreen= -1.0f, smoothRidgesBlue= -1.0f, smoothRidgesAlpha= -1.0f;

    boolean smoothValleysDrawFaded;
    float smoothValleysThreshold= -1.0f, smoothValleysRed= -1.0f, smoothValleysGreen= -1.0f, smoothValleysBlue= -1.0f, smoothValleysAlpha= -1.0f;


    boolean smoothRidgesHiddenDrawFaded;
    float smoothRidgesHiddenThreshold= -1.0f, smoothRidgesHiddenRed= -1.0f, smoothRidgesHiddenGreen= -1.0f, smoothRidgesHiddenBlue= -1.0f, smoothRidgesHiddenAlpha= -1.0f;

    boolean smoothValleysHiddenDrawFaded;    
    float smoothValleysHiddenThreshold= -1.0f, smoothValleysHiddenRed= -1.0f, smoothValleysHiddenGreen= -1.0f, smoothValleysHiddenBlue= -1.0f, smoothValleysHiddenAlpha= -1.0f;

    Vector neighbors[], adjacentFaces[];
    Vector <SegmentAndFace[][]> acrossEdge;
    Vector <float[][]> cornerareas;
    PVector pdir1[], pdir2[];
    float pointareas[], curv1[], curv2[];
    float[][] dcurv;
    
    enum TStripRep { TSTRIP_LENGTH, TSTRIP_TERM };
    Vector<Integer> tstrips;

    float boundingSphereRadius;

    float q1[][];
    float t1[][];
    float Dt1q1[][];
    float ndotv[];
    PVector normal[];

    PVector curvColors[];

    float featureSize= -1.0f, curvFeatureSize= -1.0f;

    boolean isTriangleMesh= true;

    PVector newNormal;

    /** 
     * Creates an OBJModel object and initializes it
     */
    public OBJModel() {

	renderer= (NPR)parent.g;	
	
	vertices = new Vector();
	texturev = new Vector();
	normv = new Vector();

	groups = new Hashtable();
	modelSegments = new Vector();
	materials = new Hashtable();

	debug = new Debug(parent);
	debug.enabled = true;

	minBounds= new PVector(); 
	maxBounds= new PVector();

	newNormal= new PVector(); 
    }
    
    public void dispose() {

	gl.glDeleteLists(displayList, 1);
	gl.glDeleteLists(boxyRidges, 1);
	gl.glDeleteLists(boxyValleys, 1);
	gl.glDeleteLists(boxyRidgesHidden, 1);
	gl.glDeleteLists(boxyValleysHidden, 1);
	gl.glDeleteLists(smoothRidges, 1);
	gl.glDeleteLists(smoothValleys, 1);
	gl.glDeleteLists(smoothRidgesHidden, 1);
	gl.glDeleteLists(smoothValleysHidden, 1);
	gl.glDeleteLists(apparentRidges, 1);
    }

    public PVector getMinBounds(){
	return minBounds;
    }

    public PVector getMaxBounds(){
	return maxBounds;
    }

    public float getBoundingSphereRadius(){
	return boundingSphereRadius;	
    }

    /**
     * Sets the parent applet. This is called by the renderer automatically.
     */ 
    public static void setParent(PApplet applet){

	parent= applet;
    }

    /** 
     * draws a model
     */
    public void draw() {

	if(renderer.npr()){
	    renderer.registerModel(this);
	}
	else{
	    drawModel();
	}
    }

    /**
     * Shows info on the model (number of vertices, texture coordinates, normals, groups, and segments).
     */
    public void showModelInfo() {

	debug.println("Obj Name: " + objName);
	debug.println("");
	debug.println("\tV  Size: " + vertices.size());
	debug.println("\tVt Size: " + texturev.size());
	debug.println("\tVn Size: " + normv.size());
	debug.println("");
	debug.println("\tG  Size: " + groups.size());
	debug.println("\tS Size: " + modelSegments.size());

    }


    /**
     * Disables texturing. Works only when not in NPR mode (that is, between a begin()
     * and end() call of some algorithm).
     */
    public void disableTexture() {
	flagTexture = false;
    }

    /**
     * Enables texturing. Works only when not in NPR mode (that is, between a begin()
     * and end() call of some algorithm).
     */
    public void enableTexture() {
	flagTexture = true;
    }

    // added on May 29 2007 (suggestion by polymonkey)
    /**
     * Attaches a texture to the model. Works only when not in NPR mode (that is, between a begin()
     * and end() call of some algorithm).
     * 
     * @param tex the new texture.
     */
    public void setTexture(PImage tex) {
	/*
	 * 
	 * try { PImage image = (PImage) tex.clone(); } catch (Exception e) { }
	 * if (image.width > image.height){ }
	 */
	texture = tex;

    }

    /** 
     * draws the model with Processing when not in NPR mode
     */
   
    protected void drawModel() {
	try {
	    PVector v = null, vt = null, vn = null;
	    int vtidx = 0, vnidx = 0, vidx = 0;
	    Material mtl = null;
	    // render all triangles

	    for (int s = 0; s < modelSegments.size(); s++) {

		boolean bTexture = true;

		ModelSegment tmpModelSegment = modelSegments.elementAt(s);

		mtl = (Material) materials.get(tmpModelSegment.mtlName);
		if (mtl == null) // if the material is not assigned for some
		    // reason, it uses the default material
		    // setting
		    {
			mtl = (Material) materials.get(defaultMaterialName);
			debug.println("Material '" + tmpModelSegment.mtlName
				      + "' not defined");
		    }
		if (tmpModelSegment.mtlName!=defaultMaterialName){ // this makes fill color from Processing being used
		    parent.fill(255.0f * mtl.Ka[0], 255.0f * mtl.Ka[1],
				255.0f * mtl.Ka[2], 255.0f * mtl.d);
		}
		for (int f = 0; f < tmpModelSegment.elements.size(); f++) {

		    ModelElement tmpf = (ModelElement) (tmpModelSegment.elements
							.elementAt(f));
		    parent.textureMode(PApplet.NORMALIZED);
		    parent.beginShape(mode); // specify render mode

		    if (flagTexture == false)
			bTexture = false;

		    if (tmpf.tindices.size() == 0)
			bTexture = false;

		    if(texture==null && mtl==null)
			bTexture = false;

		    if(texture==null && mtl != null && mtl.map_Kd == null) 
			bTexture = false;

		    if (bTexture)
			if (texture != null)
			    parent.texture(texture); // setting applied texture
			else
			    parent.texture(mtl.map_Kd); // setting texture from mtl info

		    if (tmpf.indices.size() > 0) {
			for (int fp = 0; fp < tmpf.indices.size(); fp++) {
			    vidx = tmpf.indices.elementAt(fp).intValue();

			    v = vertices.elementAt(vidx - 1);
			    if (v != null) {
				try {
				    if (tmpf.nindices.size() > 0) {
					vnidx = tmpf.nindices.elementAt(fp).intValue();
					vn = normv
					    .elementAt(vnidx - 1);
					parent.normal(vn.x, vn.y, vn.z);
				    }

				    if (bTexture) {
					vtidx = tmpf.tindices.elementAt(fp).intValue();

					vt = texturev
					    .elementAt(vtidx - 1);
					parent.vertex(v.x, v.y, v.z,
						      vt.x, 1.0f - vt.y);
				    } else
					parent.vertex(v.x, v.y, v.z);
				} catch (Exception e) {
				    e.printStackTrace();
				}
			    } 
			}
		    }
		    parent.endShape();
		}
	    }

	} catch (Exception e) {
	    e.printStackTrace();
	}
    }

    /**
     * Tests whether the object is a triangle mesh
     */
    protected boolean isTriangleMesh(){
	    
	for (int s = 0; s < modelSegments.size(); s++) {

	    ModelSegment tmpModelSegment = (ModelSegment) modelSegments.elementAt(s);

	    for (int f = 0; f < tmpModelSegment.elements.size(); f++) {

		ModelElement tmpf = (ModelElement) (tmpModelSegment.elements.elementAt(f));
					
		if (tmpf.indices.size() != 3) {
		    System.out.println("Not a triangle Mesh");
		    return false;
		}
	    }
	}
	return true;				  
    }

    /**
     * calculate normals
     * uses average of per-face normals, weighted according to: Max, N.
     * "Weights for Computing Vertex Normals from Facet Normals,"
     * Journal of Graphics Tools, Vol. 4, No. 2, 1999.
     */
    protected void calculateNormals() {
        
	normv= new Vector(vertices.size()); 
	for(int i= 0; i<vertices.size(); i++){
	    normv.add(new PVector());
	}

	renderer= (NPR)parent.g;
	gl= renderer.gl;

	try {
	    PVector v = null;
	    int vidx = 0, nidx= 0;

		       
	    for (int s = 0; s < modelSegments.size(); s++) {

		ModelSegment tmpModelSegment = (ModelSegment) modelSegments
		    .elementAt(s);

		for (int f = 0; f < tmpModelSegment.elements.size(); f++) {

		    ModelElement tmpf = (ModelElement) (tmpModelSegment.elements.elementAt(f));
					
		    if (tmpf.indices.size() > 0) {
						   
			PVector[] faceVertices= new PVector[3];
						
			for (int fp = 0; fp < 3; fp++) {
			    vidx= tmpf.indices.elementAt(fp).intValue();
			    tmpf.nindices.add(Integer.valueOf(vidx));
			    faceVertices[fp]= vertices.elementAt(vidx - 1);  
			}
						
			PVector a = PVector.sub(faceVertices[0], faceVertices[1]);
			PVector b = PVector.sub(faceVertices[1], faceVertices[2]);
			PVector c = PVector.sub(faceVertices[2], faceVertices[0]);
						
			PVector faceNormal= PVector.cross(a, b, null);

			float l2a = parent.sq(a.mag());
			float l2b = parent.sq(b.mag());
			float l2c = parent.sq(c.mag());
						
			nidx = tmpf.nindices.elementAt(0).intValue();
			normv.elementAt(nidx - 1).add(PVector.mult(faceNormal, 1.0f / (l2a * l2c)));

			nidx = tmpf.nindices.elementAt(1).intValue();
			normv.elementAt(nidx - 1).add(PVector.mult(faceNormal, 1.0f / (l2b * l2a)));
					
			nidx = tmpf.indices.elementAt(2).intValue();
			normv.elementAt(nidx - 1).add(PVector.mult(faceNormal, 1.0f / (l2c * l2b)));
		    }
		}
	    }


	for(int i= 0; i<vertices.size(); i++){
	    normv.elementAt(i).normalize();
	    normv.elementAt(i).mult(-1.0f);
	}

	}catch (Exception e) {
	    e.printStackTrace();
	}
    }

    /**
     * Draws the model using OpenGL commands. Used when rendering to display list.
     */
    protected void drawModelGL() {

	renderer= (NPR)parent.g;
	gl= renderer.gl;

	//if(isTriangleMesh) computeCurvColors();

	try {
	    PVector v = null, vt = null, vn = null;
	    int vtidx = 0, vnidx = 0, vidx = 0;
	    Material mtl = null;
	    // render all triangles

	    for (int s = 0; s < modelSegments.size(); s++) {

		boolean bTexture = true;

		ModelSegment tmpModelSegment = (ModelSegment) modelSegments
		    .elementAt(s);

		mtl = (Material) materials.get(tmpModelSegment.mtlName);
		if (mtl == null) // if the material is not assigned for some
		    // reason, it uses the default material
		    // setting
		    {
			mtl = (Material) materials.get(defaultMaterialName);
			debug.println("Material '" + tmpModelSegment.mtlName
				      + "' not defined");
		    }
		if (tmpModelSegment.mtlName!=defaultMaterialName){ // this makes fill color from Processing being used
		    gl.glColor4f(mtl.Ka[0], mtl.Ka[1], mtl.Ka[2], mtl.d);
		}
		for (int f = 0; f < tmpModelSegment.elements.size(); f++) {

		    ModelElement tmpf = (ModelElement) (tmpModelSegment.elements.elementAt(f));

		    // initialize bounding box 
		    if(s==0 && f==0){
			vidx = tmpf.indices.elementAt(0).intValue();
			v = vertices.elementAt(vidx - 1);
			if (v != null) {
			    minBounds.x= maxBounds.x= v.x;
			    minBounds.y= maxBounds.y= v.y;
			    minBounds.z= maxBounds.z= v.z;
			}	
		    }
		    
		    parent.textureMode(PApplet.NORMALIZED);
		    switch(mode){// specify render mode
		    case PApplet.POINTS:
			gl.glBegin(GL.GL_POINTS);
			break;
		    case PApplet.LINES:
			gl.glBegin(GL.GL_LINES);
			break;
		    case PApplet.TRIANGLES:
			gl.glBegin(GL.GL_TRIANGLES);
			break;
		    case PApplet.TRIANGLE_FAN:
			gl.glBegin(GL.GL_TRIANGLE_FAN);
			break;
		    case PApplet.TRIANGLE_STRIP:
			gl.glBegin(GL.GL_TRIANGLE_STRIP);
			break;
		    case PApplet.QUADS:
			gl.glBegin(GL.GL_QUADS);
			break;
		    case PApplet.QUAD_STRIP:
			gl.glBegin(GL.GL_QUAD_STRIP);
			break;
		    case PApplet.POLYGON:
			gl.glBegin(GL.GL_POLYGON);
			break;
		    }

		    if (flagTexture == false)
			bTexture = false;

		    if (tmpf.tindices.size() == 0)
			bTexture = false;

		    if(texture==null && mtl==null)
			bTexture = false;

		    if(texture==null && mtl != null && mtl.map_Kd == null) 
			bTexture = false;
		    /*
		      if (bTexture)
		      if (texture != null)
		      renderer.bindTexture(texture); // setting applied texture
		      else
		      renderer.bindTexture(mtl.map_Kd); // setting texture from mtl info
		    */
		    if (tmpf.indices.size() > 0) {
			for (int fp = 0; fp < tmpf.indices.size(); fp++) {
			    vidx = tmpf.indices.elementAt(fp).intValue();
			    v = vertices.elementAt(vidx - 1);
			    if (v != null) {
				try {
				    if (tmpf.nindices.size() > 0) {
					vnidx = tmpf.nindices.elementAt(fp).intValue();
					vn = normv.elementAt(vnidx - 1);
					gl.glNormal3f(vn.x, vn.y, vn.z);
				    }
				    if (tmpf.tindices.size() > 0) {
					vtidx = tmpf.tindices.elementAt(fp).intValue();

					vt = texturev.elementAt(vtidx - 1);
					gl.glTexCoord2f (vt.x, 1.0f - vt.y);
				    }
// 				    if(isTriangleMesh){
//  				    gl.glColor3f(curvColors[vidx - 1].x,
// 				    	 curvColors[vidx - 1].y,
// 					 curvColors[vidx - 1].z);
//  				    }
				    gl.glVertex3f(v.x, v.y, v.z);

				    // calculate axis aligned bounding box
				    if(v.x<minBounds.x) minBounds.x= v.x;
				    if(v.x>maxBounds.x) maxBounds.x= v.x;
				    if(v.y<minBounds.y) minBounds.y= v.y;
				    if(v.y>maxBounds.y) maxBounds.y= v.y;
				    if(v.z<minBounds.z) minBounds.z= v.z;
				    if(v.z>maxBounds.z) maxBounds.z= v.z;
				} catch (Exception e) {
				    e.printStackTrace();
				}
			    }
			}
		    }
		    gl.glEnd();
		}
	    }

	} catch (Exception e) {
	    e.printStackTrace();
	}
	boundingSphereRadius= Math.max(Math.abs(minBounds.x), 
				       Math.max(Math.abs(minBounds.y), 
						Math.max(Math.abs(minBounds.z),
							 Math.max(Math.abs(maxBounds.x), 
								  Math.max(Math.abs(maxBounds.y),
									   Math.abs(maxBounds.z))))));

    }

    /**
     * Sets the draw mode
     * 
     * @param mode one of POINTS, LINES, TRIANGLES, TRIANGLE_FAN, TRIANGLE_STRIP, QUADS, QUAD_STRIP, POLYGON.
     */
    public void drawMode(int mode) {
	this.mode = mode;
    }

    protected BufferedReader getBufferedReader(String filename) {

	BufferedReader retval = null;

	try {

	    URL url = null;
	    InputStream is = null;

	    /*
	     * parent.openStream(arg0); if (filename.startsWith("http://")) {
	     * try { url = new URL(filename); retval = new BufferedReader(new
	     * InputStreamReader(parent.openStream(filename))); return retval; }
	     * catch (MalformedURLException e) { e.printStackTrace(); return
	     * null; } catch (IOException ioe) { ioe.printStackTrace(); return
	     * null; } }
	     */

	    is = parent.openStream(filename);
	    if (is != null) {
		try {
		    retval = new BufferedReader(new InputStreamReader(is));
		    return retval;
		} catch (Exception ioe) {
		    ioe.printStackTrace();
		    return null;
		}
	    }

	    /*
	     * is = getClass().getResourceAsStream("/data/" + filename); if (is !=
	     * null) { try { retval = new BufferedReader(new
	     * InputStreamReader(is)); return retval; } catch (Exception ioe) {
	     * ioe.printStackTrace(); return null; } }
	     * 
	     * url = getClass().getResource("/" + filename); if (url != null) {
	     * System.out.println(url.toString()); try { url = new
	     * URL(filename); retval = new BufferedReader(new
	     * InputStreamReader(parent.openStream())); return retval; } catch
	     * (MalformedURLException e) { e.printStackTrace(); return null; }
	     * catch (IOException ioe) { ioe.printStackTrace(); return null; } }
	     * 
	     * url = getClass().getResource("/data/" + filename); if (url !=
	     * null) { System.out.println(url.toString()); try { url = new
	     * URL(filename); retval = new BufferedReader(new
	     * InputStreamReader(url .openStream())); return retval; } catch
	     * (MalformedURLException e) { e.printStackTrace(); return null; }
	     * catch (IOException ioe) { ioe.printStackTrace(); return null; } }
	     * 
	     * try { // look inside the sketch folder (if set) String location =
	     * parent.sketchPath + File.separator + "data"; File file = new
	     * File(location, filename); if (file.exists()) { retval = new
	     * BufferedReader(new FileReader(file)); return retval; } } catch
	     * (IOException e) { e.printStackTrace(); return null; } // ignored
	     * 
	     * try { File file = new File("data", filename); if (file.exists()) {
	     * retval = new BufferedReader(new FileReader(file)); return retval; } }
	     * catch (IOException ioe) { ioe.printStackTrace(); }
	     * 
	     * try { File file = new File(filename); if (file.exists()) { retval =
	     * new BufferedReader(new FileReader(file)); return retval; } }
	     * catch (IOException ioe) { ioe.printStackTrace(); return null; }
	     */

	} catch (SecurityException se) {
	} // online, whups

	parent.die("Could not find .OBJ file " + filename, null);

	return retval;
    }

    /**
     * load a model from a Wavefront .obj file
     *
     * @param filename the file name
     */
    public void load(String filename) {
		
	boxyRidges= -1;
	boxyValleys= -1;

	boxyRidgesThreshold= -1.0f;
	boxyRidgesRed= -1.0f;
	boxyRidgesGreen= -1.0f;
	boxyRidgesBlue= -1.0f;
	boxyRidgesAlpha= -1.0f;

	boxyValleysThreshold= -1.0f;
	boxyValleysRed= -1.0f;
	boxyValleysGreen= -1.0f;
	boxyValleysBlue= -1.0f;
	boxyValleysAlpha= -1.0f;

	smoothRidges= -1;
	smoothValleys= -1;

	smoothRidgesThreshold= -1.0f;
	smoothRidgesRed= -1.0f;
	smoothRidgesGreen= -1.0f;
	smoothRidgesBlue= -1.0f;
	smoothRidgesAlpha= -1.0f;

	smoothValleysThreshold= -1.0f;
	smoothValleysRed= -1.0f;
	smoothValleysGreen= -1.0f;
	smoothValleysBlue= -1.0f;
	smoothValleysAlpha= -1.0f;


	featureSize= -1.0f;
	curvFeatureSize= -1.0f;

	renderer= (NPR)parent.g;
	gl= renderer.gl;

	parseOBJ(getBufferedReader(filename));
	isTriangleMesh= isTriangleMesh();
	if(normv.isEmpty() && isTriangleMesh){
	    calculateNormals();
	}
	drawModelGL();	
	displayList= gl.glGenLists(1);
	gl.glNewList(displayList, GL.GL_COMPILE);
	drawModelGL();
	gl.glEndList();
	this.showModelInfo();		
    }

    /**
     * parses a Wavefront .obj file.
     */
    protected void parseOBJ(BufferedReader bread) {
	try {

	    // adding default variables to the global data table
	    // creating the default group

	    groups.put("default", defaultGroup);
	    Group currentGroup = defaultGroup;

	    // creating the default material

	    Material defaultMaterial = new Material();
	    defaultMaterial.mtlName = defaultMaterialName;
	    materials.put(defaultMaterialName, defaultMaterial);
	    String currentMaterial = defaultMaterialName;

	    // creating the default model segment

	    modelSegments.add(defaultModelSegment);
	    defaultModelSegment.mtlName = currentMaterial;
	    currentGroup.segments.add(defaultModelSegment);
	    ModelSegment currentModelSegment = defaultModelSegment;

	    String line;
	    while ((line = bread.readLine()) != null) {
		// debug.println(line);
		// parse the line
		String[] elements = line.split("\\s+");

		// if not a blank line, process the line.
		if (elements.length > 0) {

		    // analyze the format
		    if (elements[0].equals("v")) { // point vector
			PVector tmpv = new PVector();
			tmpv.x = Float.valueOf(elements[1]).floatValue();
			tmpv.y = -Float.valueOf(elements[2]).floatValue();
			tmpv.z = Float.valueOf(elements[3]).floatValue();
			vertices.add(tmpv);
		    } else if (elements[0].equals("vn")) { // normal vector
			PVector tmpv = new PVector();
			tmpv.x = Float.valueOf(elements[1]).floatValue();
			tmpv.y = -Float.valueOf(elements[2]).floatValue();
			tmpv.z = Float.valueOf(elements[3]).floatValue();
			tmpv.normalize();
			normv.add(tmpv);
		    } else if (elements[0].equals("vt")) {
			PVector tmpv = new PVector();
			tmpv.x = Float.valueOf(elements[1]).floatValue();
			tmpv.y = Float.valueOf(elements[2]).floatValue();
			texturev.add(tmpv);
		    } else if (elements[0].equals("o")) {
			if (elements[1] != null)
			    objName = elements[1];
		    } else if (elements[0].equals("mtllib")) {
			if (elements[1] != null)
			    this.parseMTL(this.getBufferedReader(elements[1]));
		    }

		    // elements that needs to consider the current context

		    else if (elements[0].equals("g")) { // grouping
			// setting
			ModelSegment newModelSegment = new ModelSegment();
			modelSegments.add(newModelSegment);
			currentModelSegment = newModelSegment;
			currentModelSegment.mtlName = currentMaterial;
						
			String groupName = "";
			for (int e = 1; e < elements.length; e++)
			    groupName += (" "+elements[e]);
						
			if (groups.get(groupName) == null) {
			    debug.println("group '" + groupName +"'");
			    Group newGroup = new Group(groupName);
			    groups.put(groupName, newGroup);
			} else {
			    debug.println("double declaration of a same group name.");
			}
						
		    } else if (elements[0].equals("usemtl")) { // material
			// setting

			String mtlName = "";
			for (int e = 1; e < elements.length; e++)
			    mtlName += (elements[e]);

			ModelSegment newModelSegment = new ModelSegment();
			modelSegments.add(newModelSegment);
			currentModelSegment = newModelSegment;
			currentModelSegment.mtlName = mtlName;
			currentMaterial = mtlName;
						
		    } else if (elements[0].equals("f")) { // Element
			ModelElement tmpf = new ModelElement();
			if (elements.length < 3) {
			    debug
				.println("Warning: potential model data error");
			}
			for (int i = 1; i < elements.length; i++) {
			    String seg = elements[i];
			    if (seg.indexOf("/") > 0) {
				String[] forder = seg.split("/");

				if (forder.length > 2) {
				    if (forder[2].length() > 0)
					tmpf.nindices.add(Integer
							  .valueOf(forder[2]));
				    if (forder[1].length() > 0)
					tmpf.tindices.add(Integer
							  .valueOf(forder[1]));
				    if (forder[0].length() > 0)
					tmpf.indices.add(Integer
							 .valueOf(forder[0]));
				} else if (forder.length > 1) {
				    if (forder[1].length() > 0)
					tmpf.tindices.add(Integer
							  .valueOf(forder[1]));
				    if (forder[0].length() > 0)
					tmpf.indices.add(Integer
							 .valueOf(forder[0]));
				} else if (forder.length > 0) {
				    if (forder[0].length() > 0)
					tmpf.indices.add(Integer
							 .valueOf(forder[0]));
				}
			    } else {
				if (seg.length() > 0)
				    tmpf.indices.add(Integer.valueOf(seg));
			    }
			}
			currentModelSegment.elements.add(tmpf);
		    } else if (elements[0].equals("ll")) { // line
			ModelElement tmpf = new ModelElement();
			tmpf.iType = ModelElement.POLYGON;
			if (elements.length < 2) {
			    debug
				.println("Warning: potential model data error");
			}
			for (int i = 1; i < elements.length; i++) {
			    tmpf.indices.add(Integer.valueOf(elements[i]));
			}
			currentModelSegment.elements.add(tmpf);
		    }
		}
	    }
	} catch (Exception e) {
	    e.printStackTrace();
	}
    }

    /**
     * Parses a material.
     */
    protected void parseMTL(BufferedReader bread) {
	try {
	    String line;
	    Material currentMtl = null;
	    while ((line = bread.readLine()) != null) {
		// parse the line
		String elements[] = line.split("\\s+");
		if (elements.length > 0) {
		    // analyze the format

		    if (elements[0].equals("newmtl")) {
			debug.println("material '" + elements[1] + "'");
			String mtlName = elements[1];
			Material tmpMtl = new Material();
			currentMtl = tmpMtl;
			materials.put(mtlName, tmpMtl);
		    } else if (elements[0].equals("map_Ka")
			       && elements.length > 1) {
			debug
			    .println("\ttexture ambient '" + elements[1]
				     + "'");
			String texname = elements[1];
			// currentMtl.map_Ka = parent.loadImage(texname);
		    } else if (elements[0].equals("map_Kd")
			       && elements.length > 1) {
			debug
			    .println("\ttexture diffuse '" + elements[1]
				     + "'");
			String texname = elements[1];
			currentMtl.map_Kd = parent.loadImage(texname);
		    } else if (elements[0].equals("Ka") && elements.length > 1) {
			currentMtl.Ka[0] = Float.valueOf(elements[1])
			    .floatValue();
			currentMtl.Ka[1] = Float.valueOf(elements[2])
			    .floatValue();
			currentMtl.Ka[2] = Float.valueOf(elements[3])
			    .floatValue();
		    } else if (elements[0].equals("d") && elements.length > 1) {
			currentMtl.d = Float.valueOf(elements[1]).floatValue();
		    }
		}
	    }
	} catch (Exception e) {
	    e.printStackTrace();
	}
    }

    /**
      * Gets the group size
      *
      * @return the group size
      */
    public int getGroupsize() {
	return this.groups.size();
    }
 
    protected Group getGroup(String groupName) {
	return (Group) this.groups.get(groupName);
    }

    /** 
     * Gets the number of vertices.
     *
     * @return the number of vertices.
     */
    public int getVertexSize() {
	return this.vertices.size();
    }

    /**
     * Return vertex with specified index.
     * 
     * @param i the index of the vertex.
     * @return the vertex.
     */
    public PVector getVertex(int i) {
	return vertices.elementAt(i);
    }

    /**
     * Sets vertex with specified index to specified vertex.
     * @param i the index of the vertex that should be set.
     * @param vertex the new values for the vertex
     */
    public void setVertex(int i, PVector vertex) {
	PVector tmpv = vertices.elementAt(i);
	tmpv.x = vertex.x;
	tmpv.y = vertex.y;
	tmpv.z = vertex.z;
    }

    /**
      * Updates the display list. Call this when done with animation.
      */
    public void update(){
        gl.glNewList(displayList, GL.GL_COMPILE);
	drawModelGL();
	gl.glEndList();
    }

    /**
     * Enable debug mode.
     */
    public void debugMode() {
	debug.enabled = true;
    }

    /**
     * Squared distance between two vectors
     */
    protected float dist2(PVector v1, PVector v2){
    
	return parent.sq(v2.x-v1.x) + 
	    parent.sq(v2.y-v1.y) + 
	    parent.sq(v2.z-v1.z);
    }


    /** 
     * Compute a "feature size" for the mesh: computed as 1% of the reciprocal of the 10-th percentile curvature
     */
    protected void needCurvFeatureSize(){

	if(!isTriangleMesh) return;
	//if(curvFeatureSize<0){
	    needCurvatures();
	    int nv = curv1.length;
	    int nsamp = Math.min(nv, 500);
	    
	    Vector<Float> samples= new Vector<Float>(nsamp * 2);
	    
	    int randq = 0;
	    for (int i = 0; i < nsamp; i++) {
		// Quick 'n dirty portable random number generator
		randq = 1664525 * randq + 1013904223;
		
		int ind = randq % nv;
		if(ind<0) ind += nv;
		samples.add(Float.valueOf(Math.abs(curv1[ind])));
		samples.add(Float.valueOf(Math.abs(curv2[ind])));
	    }
	    
	    float frac = 0.1f;
	    float mult = 0.01f;
	    float maxFeatureSize = 0.05f * boundingSphereRadius;
	    int which = (int)(frac * samples.size());
	    Collections.sort(samples);
	    curvFeatureSize= Math.min(mult / samples.elementAt(which).floatValue(), maxFeatureSize);
	    //}
    }

    /**
     * A characteristic "feature size" for the mesh.  Computed as an approximation to the median edge length
     */
    protected void needFeatureSize(){
	
	if(featureSize<0){

	    Vector<SegmentAndFace> faces= new Vector<SegmentAndFace>();
	    int numFaces= 0;
	    
	    for(int s= 0; s<modelSegments.size(); s++){
		ModelSegment tmpModelSegment = (ModelSegment) modelSegments.elementAt(s);
		numFaces+= tmpModelSegment.elements.size();
		for (int f = 0; f < tmpModelSegment.elements.size(); f++) {
		    faces.add(new SegmentAndFace(s, f));
		}
	    }
	    
	    int nf= 0;
	    int nsamp = Math.min(numFaces / 2, 333);
	    
	    Vector<Float> samples= new Vector<Float>(nsamp * 3);
	    
	    int randq= 0;
	    
	    for (int i = 0; i < nsamp; i++) {
		// Quick 'n dirty portable random number generator
		randq = 1664525 * randq + 1013904223;
		
		int ind = randq % numFaces;
		
		SegmentAndFace sf= faces.get(ind);
		int s= sf.segment, f= sf.face;
		
		ModelElement tmpf = modelSegments.elementAt(s).elements.elementAt(f);
		int v0 = (tmpf.indices.elementAt(0)).intValue();
		int v1 = (tmpf.indices.elementAt(1)).intValue();
		int v2 = (tmpf.indices.elementAt(2)).intValue();
		
		PVector p0 = vertices.elementAt(v0-1);
		PVector p1 = vertices.elementAt(v1-1);
		PVector p2 = vertices.elementAt(v2-1);
		samples.add(Float.valueOf(dist2(p0,p1)));
		samples.add(Float.valueOf(dist2(p1,p2)));
		samples.add(Float.valueOf(dist2(p2,p0)));
	    }
	    Collections.sort(samples);
	    featureSize= (float)Math.sqrt(samples.elementAt(samples.size()/2));
	}
    }
	
	
     protected float sgn(float x){
	    
	    return (x < 0) ? -1 : 1;
    }


    protected PVector hsv(float h, float s, float v){
		// From FvD
		if (s <= 0.0f)
			return new PVector(v,v,v);
		h = h % (float)(2.0f * parent.PI);
		if (h < 0.0)
			h += (float)(2.0 * parent.PI);
		h /= (float)(parent.PI / 3.0);
		int i = (int)Math.floor(h);
		float f = h - i;
		float p = v * (1.0f - s);
		float q = v * (1.0f - (s*f));
		float t = v * (1.0f - (s*(1.0f-f)));
		switch(i) {
			case 0: return new PVector(v, t, p);
			case 1: return new PVector(q, v, p);
			case 2: return new PVector(p, v, t);
			case 3: return new PVector(p, q, v);
			case 4: return new PVector(t, p, v);
			default: return new PVector(v, p, q);
		}
	}


    /**
     * Color the mesh by curvatures
     */
    protected void computeCurvColors(){
	
	needCurvFeatureSize();
	float cscale = parent.sq(8.0f * curvFeatureSize);
	
	int nv = vertices.size();
	curvColors= new PVector[nv];
	
	for (int i = 0; i < nv; i++) {
	        float H = 0.5f * (curv1[i] + curv2[i]);
		float K = curv1[i] * curv2[i];
		float h = 4.0f / 3.0f * Math.abs((float)Math.atan2(H*H-K,H*H*sgn(H)));
		float s = parent.TWO_PI * (float)Math.atan((2.0f*H*H-K)*cscale);
		curvColors[i] = hsv(h,s,1.0f);
	}
    }
    
    /**
     * Draw ridges(valleys) of a mesh by testing dihedral angles
     */
         
    public int getBoxyRidges(boolean doRidge, boolean doTest, float threshold, float red, float green, float blue, float alpha){

	if(!isTriangleMesh) return -1;

	if(doRidge){
	    if(boxyRidgesThreshold==threshold && boxyRidgesRed==red && boxyRidgesGreen==green && boxyRidgesBlue==blue && boxyRidgesAlpha==alpha){
		return boxyRidges;
	    }
	    else{
		boxyRidgesThreshold= threshold;
		boxyRidgesRed= red;
		boxyRidgesGreen= green;
		boxyRidgesBlue= blue;
		boxyRidgesAlpha= alpha;
	    }
	}
	else{
	    if(boxyValleysThreshold==threshold && boxyValleysRed==red && boxyValleysGreen==green && boxyValleysBlue==blue && boxyValleysAlpha==alpha){
	    return boxyValleys;
	    }
	    else{
		boxyValleysThreshold= threshold;
		boxyValleysRed= red;
		boxyValleysGreen= green;
		boxyValleysBlue= blue;
		boxyValleysAlpha= alpha;
	    }
	}

	needAcrossEdge();

	if(doRidge){
	    boxyRidges= gl.glGenLists(1);
	    gl.glNewList(boxyRidges, GL.GL_COMPILE);
	}
	else{
	    boxyValleys= gl.glGenLists(1);
	    gl.glNewList(boxyValleys, GL.GL_COMPILE);
	}

	gl.glBegin(GL.GL_LINES);

	PVector faceVertices1[] = new PVector[3];
	PVector faceVertices2[] = new PVector[3];
	int vidx;

	for (int s = 0; s < modelSegments.size(); s++){
	    ModelSegment tmpModelSegment = (ModelSegment) modelSegments.elementAt(s);
	    for (int f = 0; f < tmpModelSegment.elements.size(); f++) {
		ModelElement tmpf = (ModelElement) (tmpModelSegment.elements.elementAt(f));
		for (int fp = 0; fp < 3; fp++) {
		    vidx = tmpf.indices.elementAt(fp).intValue();
		    faceVertices1[fp]= vertices.elementAt(vidx - 1);  
		}
		
		PVector a = PVector.sub(faceVertices1[0], faceVertices1[1]);
		PVector b = PVector.sub(faceVertices1[1], faceVertices1[2]);
		
		PVector tmpFaceNormal= a.cross(b);
		tmpFaceNormal.normalize();

		for(int i= 0; i<3; i++){
		    SegmentAndFace sf= ((SegmentAndFace[][])acrossEdge.elementAt(s))[f][i];
		    if(sf.segment!= -1 && sf.face !=- 1){
			
			ModelSegment acrossSegment= modelSegments.elementAt(sf.segment);
			ModelElement acrossFace= acrossSegment.elements.elementAt(sf.face);

			for (int fp = 0; fp < 3; fp++) {
			    vidx = (acrossFace.indices.elementAt(fp)).intValue();
			    faceVertices2[fp]= vertices.elementAt(vidx - 1);
			}
						
			a = PVector.sub(faceVertices2[0], faceVertices2[1]);
			b = PVector.sub(faceVertices2[1], faceVertices2[2]);			
			PVector acrossFaceNormal= a.cross(b);
			acrossFaceNormal.normalize();

			if(tmpFaceNormal.dot(acrossFaceNormal) < 1.0-threshold){
			    // determine if we have a ridge or a valley here
			    // c1 = centroid of triangle 1
			    // c2 = centroid of triangle 2
			    // n1 = normal of triangle 1
			    // n2 = normal of triangle 2
			    
			    // You have a ridge iff (c1 - c2) DOT (n1 - n2) < 0

			    PVector c1= PVector.div(PVector.add(faceVertices1[0], 
								 PVector.add(faceVertices1[1],
									     faceVertices1[2])), 
						     3.0f);
			    PVector c2= PVector.div(PVector.add(faceVertices2[0], 
								 PVector.add(faceVertices2[1],
									     faceVertices2[2])), 
						     3.0f);

			    boolean isRidge= PVector.sub(c1, c2).dot(PVector.sub(tmpFaceNormal, acrossFaceNormal)) < 0;

			    if((doRidge && isRidge) || (!doRidge && !isRidge)){
				for (int fp = 0; fp < 3; fp++){
				    vidx = tmpf.indices.elementAt(fp).intValue();
				    if(acrossFace.indices.indexOf(Integer.valueOf(vidx)) != -1){
					PVector vertex= vertices.elementAt(vidx-1);
					gl.glColor4f(red, green, blue, alpha);
					
					gl.glVertex3f(vertex.x, vertex.y, vertex.z);  
				    }
				}
			    }
			}
		    }
		}
	    }
	}
	
	gl.glEnd();
	gl.glEndList();

	if(doRidge) return boxyRidges;
	else return boxyValleys;
    }

public int getBoxyRidgesHidden(boolean doRidge, boolean doTest, float threshold, float red, float green, float blue, float alpha){

	if(!isTriangleMesh) return -1;

	if(doRidge){
	    if(boxyRidgesHiddenThreshold==threshold && boxyRidgesHiddenRed==red && boxyRidgesHiddenGreen==green && boxyRidgesHiddenBlue==blue && boxyRidgesHiddenAlpha==alpha){
		return boxyRidgesHidden;
	    }
	    else{
		boxyRidgesHiddenThreshold= threshold;
		boxyRidgesHiddenRed= red;
		boxyRidgesHiddenGreen= green;
		boxyRidgesHiddenBlue= blue;
		boxyRidgesHiddenAlpha= alpha;
	    }
	}
	else{
	    if(boxyValleysHiddenThreshold==threshold && boxyValleysHiddenRed==red && boxyValleysHiddenGreen==green && boxyValleysHiddenBlue==blue && boxyValleysHiddenAlpha==alpha){
	    return boxyValleysHidden;
	    }
	    else{
		boxyValleysHiddenThreshold= threshold;
		boxyValleysHiddenRed= red;
		boxyValleysHiddenGreen= green;
		boxyValleysHiddenBlue= blue;
		boxyValleysHiddenAlpha= alpha;
	    }
	}

	needAcrossEdge();

	if(doRidge){
	    boxyRidgesHidden= gl.glGenLists(1);
	    gl.glNewList(boxyRidgesHidden, GL.GL_COMPILE);
	}
	else{
	    boxyValleysHidden= gl.glGenLists(1);
	    gl.glNewList(boxyValleysHidden, GL.GL_COMPILE);
	}

	gl.glBegin(GL.GL_LINES);

	PVector faceVertices1[] = new PVector[3];
	PVector faceVertices2[] = new PVector[3];
	int vidx;

	for (int s = 0; s < modelSegments.size(); s++){
	    ModelSegment tmpModelSegment = (ModelSegment) modelSegments.elementAt(s);
	    for (int f = 0; f < tmpModelSegment.elements.size(); f++) {
		ModelElement tmpf = (ModelElement) (tmpModelSegment.elements.elementAt(f));
		for (int fp = 0; fp < 3; fp++) {
		    vidx = tmpf.indices.elementAt(fp).intValue();
		    faceVertices1[fp]= vertices.elementAt(vidx - 1);  
		}
		
		PVector a = PVector.sub(faceVertices1[0], faceVertices1[1]);
		PVector b = PVector.sub(faceVertices1[1], faceVertices1[2]);
		
		PVector tmpFaceNormal= a.cross(b);
		tmpFaceNormal.normalize();

		for(int i= 0; i<3; i++){
		    SegmentAndFace sf= ((SegmentAndFace[][])acrossEdge.elementAt(s))[f][i];
		    if(sf.segment!= -1 && sf.face !=- 1){
			
			ModelSegment acrossSegment= modelSegments.elementAt(sf.segment);
			ModelElement acrossFace= acrossSegment.elements.elementAt(sf.face);

			for (int fp = 0; fp < 3; fp++) {
			    vidx = (acrossFace.indices.elementAt(fp)).intValue();
			    faceVertices2[fp]= vertices.elementAt(vidx - 1);
			}
						
			a = PVector.sub(faceVertices2[0], faceVertices2[1]);
			b = PVector.sub(faceVertices2[1], faceVertices2[2]);			
			PVector acrossFaceNormal= a.cross(b);
			acrossFaceNormal.normalize();

			if(tmpFaceNormal.dot(acrossFaceNormal) < 1.0-threshold){
			    // determine if we have a ridge or a valley here
			    // c1 = centroid of triangle 1
			    // c2 = centroid of triangle 2
			    // n1 = normal of triangle 1
			    // n2 = normal of triangle 2
			    
			    // You have a ridge iff (c1 - c2) DOT (n1 - n2) < 0

			    PVector c1= PVector.div(PVector.add(faceVertices1[0], 
								 PVector.add(faceVertices1[1],
									     faceVertices1[2])), 
						     3.0f);
			    PVector c2= PVector.div(PVector.add(faceVertices2[0], 
								 PVector.add(faceVertices2[1],
									     faceVertices2[2])), 
						     3.0f);

			    boolean isRidge= PVector.sub(c1, c2).dot(PVector.sub(tmpFaceNormal, acrossFaceNormal)) < 0;

			    if((doRidge && isRidge) || (!doRidge && !isRidge)){
				for (int fp = 0; fp < 3; fp++){
				    vidx = tmpf.indices.elementAt(fp).intValue();
				    if(acrossFace.indices.indexOf(Integer.valueOf(vidx)) != -1){
					PVector vertex= vertices.elementAt(vidx-1);
					gl.glColor4f(red, green, blue, alpha);
					
					gl.glVertex3f(vertex.x, vertex.y, vertex.z);  
				    }
				}
			    }
			}
		    }
		}
	    }
	}
	
	gl.glEnd();
	gl.glEndList();

	if(doRidge) return boxyRidgesHidden;
	else return boxyValleysHidden;
    }
    

    /**
     * Draw ridges (valleys) of the mesh
     */
    
    public int getSmoothRidges(boolean doRidge, boolean doTest, boolean drawFaded, float threshold, float red, float green, float blue, float alpha){

	if(!isTriangleMesh) return -1;


	if(doRidge){
	    if(smoothRidgesDrawFaded==drawFaded && smoothRidgesThreshold==threshold && smoothRidgesRed==red && smoothRidgesGreen==green && smoothRidgesBlue==blue && smoothRidgesAlpha==alpha){
		return smoothRidges;
	    }
	    else{
		smoothRidgesDrawFaded= drawFaded;
		smoothRidgesThreshold= threshold;
		smoothRidgesRed= red;
		smoothRidgesGreen= green;
		smoothRidgesBlue= blue;
		smoothRidgesAlpha= alpha;
	    }
	}
	else{
	    if(smoothValleysDrawFaded==drawFaded && smoothValleysThreshold==threshold && smoothValleysRed==red && smoothValleysGreen==green && smoothValleysBlue==blue && smoothValleysAlpha==alpha){
		return smoothValleys;
	    }
	    else{
		smoothValleysDrawFaded= drawFaded;
		smoothValleysThreshold= threshold;
		smoothValleysRed= red;
		smoothValleysGreen= green;
		smoothValleysBlue= blue;
		smoothValleysAlpha= alpha;
	    }
	}

	needDCurv();
	needTStrips();

	needCurvFeatureSize();

	if(doRidge){
	    smoothRidges= gl.glGenLists(1);
	    gl.glNewList(smoothRidges, GL.GL_COMPILE);
	}
	else{
	    smoothValleys= gl.glGenLists(1);
	    gl.glNewList(smoothValleys, GL.GL_COMPILE);
	}
	
	gl.glBegin(GL.GL_LINES);

	int t = 0;
	int stripend = 0;
	int end = tstrips.size();

	// Walk through triangle strips
	while (true) {
	    if (t >= stripend) {
		if (t >= end) break;
		// New strip: each strip is stored as
		// length followed by indices
		stripend = t + 1 + tstrips.get(t).intValue();
		// Skip over length plus first two indices of
		// first face
		t+= 3;
	    }
	    drawFaceRidges(tstrips.get(t-2).intValue(), 
			   tstrips.get(t-1).intValue(), 
			   tstrips.get(t).intValue(),
			   doRidge, doTest, drawFaded, 
			   threshold/curvFeatureSize, 
			   red, green, blue, alpha);
	    t++;
	}
	
	gl.glEnd();
	gl.glEndList();
	
	if(doRidge) return smoothRidges;
	else return smoothValleys;
    } 
    

    public int getSmoothRidgesHidden(boolean doRidge, boolean doTest, boolean drawFaded, float threshold, float red, float green, float blue, float alpha){

	if(!isTriangleMesh) return -1;


	if(doRidge){
	    if(smoothRidgesHiddenDrawFaded==drawFaded && smoothRidgesHiddenThreshold==threshold && smoothRidgesHiddenRed==red && smoothRidgesHiddenGreen==green && smoothRidgesHiddenBlue==blue && smoothRidgesHiddenAlpha==alpha){
		return smoothRidgesHidden;
	    }
	    else{
		smoothRidgesHiddenDrawFaded= drawFaded;
		smoothRidgesHiddenThreshold= threshold;
		smoothRidgesHiddenRed= red;
		smoothRidgesHiddenGreen= green;
		smoothRidgesHiddenBlue= blue;
		smoothRidgesHiddenAlpha= alpha;
	    }
	}
	else{
	    if(smoothValleysHiddenDrawFaded==drawFaded && smoothValleysHiddenThreshold==threshold && smoothValleysHiddenRed==red && smoothValleysHiddenGreen==green && smoothValleysHiddenBlue==blue && smoothValleysHiddenAlpha==alpha){
	    return smoothValleysHidden;
	    }
	    else{
		smoothValleysHiddenDrawFaded= drawFaded;
		smoothValleysHiddenThreshold= threshold;
		smoothValleysHiddenRed= red;
		smoothValleysHiddenGreen= green;
		smoothValleysHiddenBlue= blue;
		smoothValleysHiddenAlpha= alpha;
	    }
	}

	needDCurv();
	needTStrips();

	needCurvFeatureSize();

	if(doRidge){
	    smoothRidgesHidden= gl.glGenLists(1);
	    gl.glNewList(smoothRidgesHidden, GL.GL_COMPILE);
	}
	else{
	    smoothValleysHidden= gl.glGenLists(1);
	    gl.glNewList(smoothValleysHidden, GL.GL_COMPILE);
	}
	
	gl.glBegin(GL.GL_LINES);

	int t = 0;
	int stripend = 0;
	int end = tstrips.size();

	// Walk through triangle strips
	while (true) {
	    if (t >= stripend) {
		if (t >= end) break;
		// New strip: each strip is stored as
		// length followed by indices
		stripend = t + 1 + tstrips.get(t).intValue();
		// Skip over length plus first two indices of
		// first face
		t+= 3;
	    }
	    drawFaceRidges(tstrips.get(t-2).intValue(), 
			   tstrips.get(t-1).intValue(), 
			   tstrips.get(t).intValue(),
			   doRidge, doTest, drawFaded, 
			   threshold/curvFeatureSize, 
			   red, green, blue, alpha);
	    t++;
	}
	
	gl.glEnd();
	gl.glEndList();
	
	if(doRidge) return smoothRidgesHidden;
	else return smoothValleysHidden;
    } 

    /** Draw ridges or valleys (depending on do_ridge) in a triangle v0,v1,v2
     * - do_test checks for curvature maxima/minina for ridges/valleys
     * (when off, it draws positive minima and negative maxima)
     */

    // Algorithm based on formulas of Ohtake et al., 2004.
    protected void drawFaceRidges(int v0, int v1, int v2,
			boolean doRidge, boolean doTest, boolean drawFaded, 
			float threshold,
			float red, float green, float blue, float alpha){

	// Check if ridge possible at vertices just based on curvatures
	if(doRidge){
	    if ((curv1[v0-1] <= 0.0f) ||
		(curv1[v1-1] <= 0.0f) ||
		(curv1[v2-1] <= 0.0f)){
		return;
	    }
	}else{
	    if ((curv1[v0-1] >= 0.0f) ||
		(curv1[v1-1] >= 0.0f) ||
		(curv1[v2-1] >= 0.0f))
		return;
	}

	// Sign of curvature on ridge/valley
	float rvSign = doRidge ? 1.0f : -1.0f;

	// The "tmax" are the principal directions of maximal curvature,
	// flipped to point in the direction in which the curvature
	// is increasing (decreasing for valleys).  Note that this
	// is a bit different from the notation in Ohtake et al.,
	// but the tests below are equivalent.
	float emax0= dcurv[v0-1][0];
	float emax1= dcurv[v1-1][0]; 
	float emax2= dcurv[v2-1][0];
	PVector tmax0= PVector.mult(pdir1[v0-1], 
				    rvSign * dcurv[v0-1][0]);
	PVector tmax1= PVector.mult(pdir1[v1-1],
				    rvSign * dcurv[v1-1][0]);
	PVector tmax2= PVector.mult(pdir1[v2-1],
				    rvSign * dcurv[v2-1][0]);

	// We have a "zero crossing" if the tmaxes along an edge
	// point in opposite directions
	boolean z01= (tmax0.dot(tmax1) <= 0.0f);
	boolean z12= (tmax1.dot(tmax2) <= 0.0f);
	boolean z20= (tmax2.dot(tmax0) <= 0.0f);

	if (!((z01 && z12) || (z01 && z20) || (z12 && z20)))
	    return;

	if (doTest) {
	    PVector p0= vertices.elementAt(v0-1),
		p1= vertices.elementAt(v1-1),
		p2= vertices.elementAt(v2-1);
		
	    // Check whether we have the correct flavor of extremum:
	    // Is the curvature increasing along the edge?
	    z01 = z01 && ((tmax0.dot(PVector.sub(p1, p0))) >= 0.0f ||
			  (tmax1.dot(PVector.sub(p1, p0))) <= 0.0f);
	    z12 = z12 && ((tmax1.dot(PVector.sub(p2, p1))) >= 0.0f ||
			  (tmax2.dot(PVector.sub(p2, p1))) <= 0.0f);
	    z20 = z20 && ((tmax2.dot(PVector.sub(p0, p2))) >= 0.0f ||
			  (tmax0.dot(PVector.sub(p0, p2))) <= 0.0f);

	    if (!((z01 && z12) || (z01 && z20) || (z12 && z20)))
		return;
	}

	// Draw line segment
	float kmax0= curv1[v0-1];
	float kmax1= curv1[v1-1];
	float kmax2= curv1[v2-1];

	if (!z01) {
	    drawSegmentRidge(v1, v2, v0,
			     emax1, emax2, emax0,
			     kmax1, kmax2, kmax0,
			     threshold, false, drawFaded,
			     red, green, blue, alpha);
	} else if (!z12) {
	    drawSegmentRidge(v2, v0, v1,
			     emax2, emax0, emax1,
			     kmax2, kmax0, kmax1,
			     threshold, false, drawFaded,
			     red, green, blue, alpha);
	} else if (!z20) {
	    drawSegmentRidge(v0, v1, v2,
			     emax0, emax1, emax2,
			     kmax0, kmax1, kmax2,
			     threshold, false, drawFaded,
			     red, green, blue, alpha);
	} else {
	    // All three edges have crossings -- connect all to center
	    drawSegmentRidge(v1, v2, v0,
			     emax1, emax2, emax0,
			     kmax1, kmax2, kmax0,
			     threshold, true, drawFaded,
			     red, green, blue, alpha);
	    drawSegmentRidge(v2, v0, v1,
			     emax2, emax0, emax1,
			     kmax2, kmax0, kmax1,
			     threshold, true, drawFaded,
			     red, green, blue, alpha);
	    drawSegmentRidge(v0, v1, v2,
			     emax0, emax1, emax2,
			     kmax0, kmax1, kmax2,
			     threshold, true, drawFaded,
			     red, green, blue, alpha);
	}
    }

    /** 
     * Draw part of a ridge/valley curve on one triangle face.  
     * v0, v1, v2 are the indices of the 3 vertices; this function assumes that the curve connects points on the edges v0-v1 and v1-v2
     * (or connects point on v0-v1 to center if to_center is true)
     */
    protected void drawSegmentRidge(int v0, int v1, int v2,
			    float emax0, float emax1, float emax2,
			    float kmax0, float kmax1, float kmax2,
			    float threshold, 
			    boolean toCenter, boolean drawFaded,
			    float red, float green, float blue, float alpha){

	// Interpolate to find ridge/valley line segment endpoints
	// in this triangle and the curvatures there
	float w10 = Math.abs(emax0) / (Math.abs(emax0) + Math.abs(emax1));
	float w01 = 1.0f - w10;
	PVector p01 = PVector.add(PVector.mult(vertices.elementAt(v0-1), w01),
				  PVector.mult(vertices.elementAt(v1-1), w10));
	float k01 = Math.abs(w01 * kmax0 + w10 * kmax1);
	
	PVector p12= new PVector();
	float k12;
	if (toCenter) {
	    // Connect first point to center of triangle
	    p12 = PVector.div(PVector.add(vertices.elementAt(v0-1),
					  PVector.add(vertices.elementAt(v1-1),
						      vertices.elementAt(v2-1))), 3.0f);
	    k12 = Math.abs(kmax0 + kmax1 + kmax2) / 3.0f;
	} 
	else {
	    // Connect first point to second one (on next edge)
	    float w21 = Math.abs(emax1) / (Math.abs(emax1) + Math.abs(emax2));
	    float w12 = 1.0f - w21;
	    p12 = PVector.add(PVector.mult(vertices.elementAt(v1-1), w12),
			      PVector.mult(vertices.elementAt(v2-1), w21));
	    k12 = Math.abs(w12 * kmax1 + w21 * kmax2);
	}
	
	// Don't draw below threshold
	k01 -= threshold;
	if (k01 < 0.0f)
	    k01 = 0.0f;
	k12 -= threshold;
	if (k12 < 0.0f)
	    k12 = 0.0f;
	    
	// Skip lines that you can't see...
	if (k01 == 0.0f && k12 == 0.0f)
	    return;
	    
	// Fade lines
	if (drawFaded) {
	    k01 /= (k01 + threshold);
	    k12 /= (k12 + threshold);
	} else {
	    k01 = k12 = 1.0f;
	}
	    
	// Draw the line segment
	gl.glColor4f(red, green, blue, k01*alpha);
	gl.glVertex3f(p01.x, p01.y, p01.z);
	gl.glColor4f(red, green, blue, k12*alpha);
	gl.glVertex3f(p12.x, p12.y, p12.z);
    }


    /*
      Original by Tilke Judd
      Tweaks by Szymon Rusinkiewicz
      
      Compute apparent ridges.
      
      Implements method of
      Judd, T., Durand, F, and Adelson, E.
      Apparent Ridges for Line Drawing,
      ACM Trans. Graphics (Proc. SIGGRAPH), vol. 26, no. 3, 2007.
    */
    
    /**
     * Compute largest eigenvalue and associated eigenvector of a symmetric 2x2 matrix.  Solves characteristic equation.
     * Inputs: three elements of matrix (upper-left, diag, lower-right)
     * Outputs: largest (in magnitude) eigenvector/value
     */
    protected void largestEig2x2(float m1, float m12, float m2, float e1[], float[] l1){
	l1[0] = 0.5f * (m1 + m2);
	// The result of the below sqrt is positive, so to get the largest
	// eigenvalue we add it if we were positive already, else subtract
	if (l1[0] > 0.0f){
	    l1[0] += (float)Math.sqrt(parent.sq(m12) + 0.25f * parent.sq(m2-m1));
	}
	else{
	    l1[0] -= (float)Math.sqrt(parent.sq(m12) + 0.25f * parent.sq(m2-m1));
	}
	// Find corresponding eigenvector
	e1[0]= m2 - l1[0];
	e1[1]= -m12;
	//normalize
	float length = (float)Math.sqrt(parent.sq(e1[0])+parent.sq(e1[1]));
	e1[0]/= length;
	e1[1]/= length;
    }


    
    /**
     * Compute principal view-dependent curvatures and directions at vertex i.
     * ndotv = cosine of angle between normal and view direction
     * (u,v) = coordinates of w (projected view) in principal coordinates
     *  Pass in u^2, u*v, and v^2, since those are readily available.
     * Fills in q1 and t1 (using the paper's notation).
     * Note that the latter is expressed in the (pdir1,pdir2) coordinate basis
     */
    
    protected void computeViewdepCurv(int i, float ndotv,
			    float u2, float uv, float v2,
			    float q1[], float t1[])
    {
	// Find the entries in Q = S * P^-1
	//                       = S + (sec theta - 1) * S * w * w^T
	float sectheta_minus1 = 1.0f / Math.abs(ndotv) - 1.0f;
	float Q11= curv1[i] * (1.0f + sectheta_minus1 * u2);
	float Q12= curv1[i] * (       sectheta_minus1 * uv);
	float Q21= curv2[i] * (       sectheta_minus1 * uv);
	float Q22= curv2[i] * (1.0f + sectheta_minus1 * v2);
	
	// Find the three entries in the (symmetric) matrix Q^T Q
	float QTQ1  = Q11 * Q11 + Q21 * Q21;
	float QTQ12 = Q11 * Q12 + Q21 * Q22;
	float QTQ2  = Q12 * Q12 + Q22 * Q22;
	
	// Compute eigenstuff
	largestEig2x2(QTQ1, QTQ12, QTQ2, t1, q1);
    }

    /**
     *  Compute D_{t_1} q_1 - the derivative of max view-dependent curvature in the principal max view-dependent curvature direction.
     */
    
    protected void computeDt1q1(int vidx, float ndotv,
		      float q1[][], float[][] t1,
		      float Dt1q1[]){
	
	needAdjacentFaces();

	PVector v0 = vertices.elementAt(vidx);
	float this_viewdep_curv = q1[vidx][0];
	PVector world_t1 = PVector.add(PVector.mult(pdir1[vidx], t1[vidx][0]),
				       PVector.mult(pdir2[vidx], t1[vidx][1]));
	PVector world_t2 = world_t1.cross(normal[vidx]);
	float v0_dot_t2 = v0.dot(world_t2);
	
	Dt1q1[0] = 0.0f;
	int n = 0;
	
	int naf = adjacentFaces[vidx].size();
	for (int j = 0; j < naf; j++) {
	    // We're in a triangle adjacent to the vertex of interest.
	    // The current vertex is v0 - let v1 and v2 be the other two
	    SegmentAndFace sf = (SegmentAndFace) adjacentFaces[vidx].elementAt(j);
	    int s= sf.segment, f= sf.face;
	    int ind = modelSegments.elementAt(s).elements.elementAt(f).indices.indexOf(Integer.valueOf(vidx+1));
	    int prevInd= ((ind)>0 ? (ind)-1 : (ind)+2);
	    int i1 = modelSegments.elementAt(s).elements.elementAt(f).indices.elementAt((ind+1)%3);
	    int i2 = modelSegments.elementAt(s).elements.elementAt(f).indices.elementAt(prevInd);
	    PVector v1 = vertices.elementAt(i1-1);
	    PVector v2 = vertices.elementAt(i2-1);

	    // Find the point p on the segment between v1 and v2 such that
	    // its vector from v0 is along t1, i.e. perpendicular to t2.
	    // Linear combination: p = w1*v1 + w2*v2, where w2 = 1-w1
	    float v1_dot_t2 = v1.dot(world_t2);
	    float v2_dot_t2 = v2.dot(world_t2);
	    float w1 = (v2_dot_t2 - v0_dot_t2) / (v2_dot_t2 - v1_dot_t2);
	    
	    // If w1 is not in [0..1) then we're not interested.  
	    // Incidentally, the computation of w1 can result in infinity,
	    // but the comparison should do the right thing...
	    if (w1 < 0.0f || w1 >= 1.0f)
		continue;
		
	    // Construct the opposite point
	    float w2 = 1.0f - w1;
	    PVector p = PVector.add(PVector.mult(v1, w1), 
				    PVector.mult(v2, w2));

	    // And interpolate to find the view-dependent curvature at that point
	    float interp_viewdep_curv = w1 * q1[i1-1][0] + w2 * q1[i2-1][0];

	    // Finally, take the *projected* view-dependent curvature derivative
	    float proj_dist = PVector.sub(p, v0).dot(world_t1);
	    proj_dist *= Math.abs(ndotv);
	    Dt1q1[0] += (interp_viewdep_curv - this_viewdep_curv) / proj_dist;
	    n++;

	    // To save time, quit as soon as we have two estimates
	    // (that's all we're going to get, anyway)
	    if (n == 2) {
		Dt1q1[0] *= 0.5f;
		return;
	    }
	}
    }
    
    protected PVector trinorm(PVector v0, PVector v1, PVector v2){
	
	return PVector.mult(PVector.sub(v2, v0).cross(PVector.sub(v1, v0)), 0.5f);
    }

    /**
     * Draw part of an apparent ridge/valley curve on one triangle face.
     * v0,v1,v2 are the indices of the 3 vertices; this function assumes that the curve connects points on the edges v0-v1 and v1-v2
     * (or connects point on v0-v1 to center if to_center is true)
     */
    protected void drawSegmentApparentRidge(int v0, int v1, int v2,
					    float emax0, float emax1, float emax2,
					    float kmax0, float kmax1, float kmax2,
					    PVector tmax0, PVector tmax1, PVector tmax2,
					    float threshold, 
					    boolean toCenter, boolean doTest, boolean drawFaded,   
					    float red, float green, float blue, float alpha){
	
	// Interpolate to find ridge/valley line segment endpoints
	// in this triangle and the curvatures there
	float w10 = Math.abs(emax0) / (Math.abs(emax0) + Math.abs(emax1));
	float w01 = 1.0f - w10;
	PVector p01 = PVector.add(PVector.mult(vertices.elementAt(v0-1), w01),
				  PVector.mult(vertices.elementAt(v1-1), w10));
	float k01 = Math.abs(w01 * kmax0 + w10 * kmax1);

	PVector p12;
	float k12;
	if (toCenter){
		// Connect first point to center of triangle
	    p12 = PVector.div(PVector.add(vertices.elementAt(v0-1), 
					  PVector.add(vertices.elementAt(v1-1), 
						      vertices.elementAt(v2-1))), 3.0f);
	    k12 = Math.abs(kmax0 + kmax1 + kmax2) / 3.0f;
	} 
	else{
	    // Connect first point to second one (on next edge)
	    float w21 = Math.abs(emax1) / (Math.abs(emax1) + Math.abs(emax2));
	    float w12 = 1.0f - w21;
	    p12 = PVector.add(PVector.mult(vertices.elementAt(v1-1), w12),
			      PVector.mult(vertices.elementAt(v2-1), w21));
	    k12 = Math.abs(w12 * kmax1 + w21 * kmax2);
	}

	// Don't draw below threshold
	k01 -= threshold;
	if (k01 < 0.0f)
	    k01 = 0.0f;
	k12 -= threshold;
	if (k12 < 0.0f)
	    k12 = 0.0f;

	// Skip lines that you can't see...
	if (k01 == 0.0f && k12 == 0.0f)
	    return;

	// Perform test: do the tmax-es point *towards* the segment? (Fig 6)
	if (doTest) {
	    // Find the vector perpendicular to the segment (p01 <-> p12)
	    PVector perp = PVector.sub(p01, p12).cross(trinorm(vertices.elementAt(v0-1),
							       vertices.elementAt(v1-1),
							       vertices.elementAt(v2-1)));
	    // We want tmax1 to point opposite to perp, and
	    // tmax0 and tmax2 to point along it.  Otherwise, exit out.
	    if (tmax0.dot(perp) <= 0.0f ||
		tmax1.dot(perp) >= 0.0f ||
		tmax2.dot(perp) <= 0.0f)
		return;
	}
	
	// Fade lines
	if (drawFaded) {
	    k01 /= (k01 + threshold);
	    k12 /= (k12 + threshold);
	}else {
	    k01 = k12 = 1.0f;
	}
	
	// Draw the line segment
       
	gl.glColor4f(red, green, blue, k01*alpha);
	gl.glVertex3f(p01.x, p01.y, p01.z);
	gl.glColor4f(red, green, blue, k12*alpha);
	gl.glVertex3f(p12.x, p12.y, p12.z);
    }
    


    /**
     * Draw apparent ridges in a triangle
     */
     protected void drawFaceApparentRidges(int v0, int v1, int v2,
					   float ndotv[], float q1[][],
					   float t1[][], float Dt1q1[][],
					   boolean doTest, boolean drawFaded, float threshold,
					   float red, float green, float blue, float alpha){

	// Trivial reject if this face isn't getting past the threshold anyway
	float kmax0 = q1[v0-1][0];
	float kmax1 = q1[v1-1][0];
	float kmax2 = q1[v2-1][0];
	if (kmax0 <= threshold && kmax1 <= threshold && kmax2 <= threshold)
	    return;
	
	// The "tmax" are the principal directions of view-dependent curvature,
	// flipped to point in the direction in which the curvature
	// is increasing.
	float emax0 = Dt1q1[v0-1][0];
	float emax1 = Dt1q1[v1-1][0];
	float emax2 = Dt1q1[v2-1][0];
	PVector world_t1_0 = PVector.add(PVector.mult(pdir1[v0-1], t1[v0-1][0]), 
					 PVector.mult(pdir2[v0-1], t1[v0-1][1]));
	PVector world_t1_1 = PVector.add(PVector.mult(pdir1[v1-1], t1[v1-1][0]),
					 PVector.mult(pdir2[v1-1], t1[v1-1][1]));
	PVector world_t1_2 =  PVector.add(PVector.mult(pdir1[v2-1], t1[v2-1][0]),
					  PVector.mult(pdir2[v2-1], t1[v2-1][1]));
	PVector  tmax0 = PVector.mult(world_t1_0, Dt1q1[v0-1][0]);
	PVector  tmax1 = PVector.mult(world_t1_1, Dt1q1[v1-1][0]);
	PVector  tmax2 = PVector.mult(world_t1_2, Dt1q1[v2-1][0]);

	// We have a "zero crossing" if the tmaxes along an edge
	// point in opposite directions
	boolean z01 = (tmax0.dot(tmax1) <= 0.0f);
	boolean z12 = (tmax1.dot(tmax2) <= 0.0f);
	boolean z20 = (tmax2.dot(tmax0) <= 0.0f);

	if (!((z01 && z12) || (z01 && z20) || (z12 && z20)))
	    return;

	// Draw line segment
	if (!z01) {
	    drawSegmentApparentRidge(v1, v2, v0,
				     emax1, emax2, emax0,
				     kmax1, kmax2, kmax0,
				     tmax1, tmax2, tmax0,
				     threshold, false, doTest, drawFaded,
				     red, green, blue, alpha);
	} else if (!z12) {
	    drawSegmentApparentRidge(v2, v0, v1,
				     emax2, emax0, emax1,
				     kmax2, kmax0, kmax1,
				     tmax2, tmax0, tmax1,
				     threshold, false, doTest, drawFaded,
				     red, green, blue, alpha);
	} else if (!z20) {
	    drawSegmentApparentRidge(v0, v1, v2,
				     emax0, emax1, emax2,
				     kmax0, kmax1, kmax2,
				     tmax0, tmax1, tmax2,
				     threshold, false, doTest, drawFaded,
				     red, green, blue, alpha);
	} else {
	    // All three edges have crossings -- connect all to center
	    drawSegmentApparentRidge(v1, v2, v0,
				     emax1, emax2, emax0,
				     kmax1, kmax2, kmax0,
				     tmax1, tmax2, tmax0,
				     threshold, true, doTest, drawFaded,
				     red, green, blue, alpha);
	    drawSegmentApparentRidge(v2, v0, v1,
				     emax2, emax0, emax1,
				     kmax2, kmax0, kmax1,
				     tmax2, tmax0, tmax1,
				     threshold, true, doTest, drawFaded,
				     red, green, blue, alpha);
	    drawSegmentApparentRidge(v0, v1, v2,
				     emax0, emax1, emax2,
				     kmax0, kmax1, kmax2,
				     tmax0, tmax1, tmax2,
				     threshold, true, doTest, drawFaded,
				     red, green, blue, alpha);
	}
     }


    /* Draw apparent ridges of the mesh */
    protected void drawApparentRidges(PVector viewPos, 
				      boolean doTest, boolean drawFaded, float threshold,
				      float red, float green, float blue, float alpha){
	
	if(!isTriangleMesh) return;

	needCurvFeatureSize();
	needTStrips();

	int nv= vertices.size();
	if(ndotv== null || ndotv.length != nv){
	    ndotv= new float[nv];
	}
	if(q1== null || q1.length!=nv){

	    q1= new float[nv][1];
	    t1= new float[nv][2];
	    Dt1q1= new float[nv][1];
	    if(normal==null || normal.length!=vertices.size()){
		normal= new PVector[nv];
		for (int s = 0; s < modelSegments.size(); s++){
		    ModelSegment tmpModelSegment = (ModelSegment) modelSegments.elementAt(s);
		    for (int f = 0; f < tmpModelSegment.elements.size(); f++) {
			ModelElement tmpf = (ModelElement) (tmpModelSegment.elements.elementAt(f));
			for(int i= 0; i<3; i++){
			    int v= tmpf.indices.elementAt(i).intValue();
			    int n= tmpf.nindices.elementAt(i).intValue();
			    normal[v-1]= normv.elementAt(n-1);
			}
		    }
		}
	    }    
	}

	for(int i= 0; i<nv; i++){
				  
	    PVector viewDir = PVector.sub(viewPos, vertices.elementAt(i));
	    float rlv = 1.0f / viewDir.mag();
	    viewDir.mult(rlv);
	    ndotv[i]= viewDir.dot(normal[i]);

	    float u = viewDir.dot(pdir1[i]), u2 = u*u;
	    float v = viewDir.dot(pdir2[i]), v2 = v*v;

	    float csc2theta = 1.0f / (u2 + v2);
	    computeViewdepCurv(i, ndotv[i],
				 u2*csc2theta, u*v*csc2theta, v2*csc2theta,
				 q1[i], t1[i]);
	}
	
	for (int i = 0; i < nv; i++)
	    computeDt1q1(i, ndotv[i], q1, t1, Dt1q1[i]);

	gl.glBegin(GL.GL_LINES);

	int t = 0;
	int stripend = 0;
	int end = tstrips.size();

	// Walk through triangle strips
	while (true) {
	    if (t >= stripend) {
		if (t >= end)
		    break;
		// New strip: each strip is stored as
		// length followed by indices
		stripend = t + 1 + tstrips.get(t);
		// Skip over length plus first two indices of
		// first face
		t += 3;
	    }
	    drawFaceApparentRidges(tstrips.get(t-2).intValue(),
				   tstrips.get(t-1).intValue(),
				   tstrips.get(t).intValue(),
				   ndotv, q1, t1, Dt1q1,
				   doTest, drawFaded, (threshold/parent.sq(curvFeatureSize)),
				   red, green, blue, alpha);
	    t++;
	}

	gl.glEnd();
    }

    /* end apparent ridges */



    /**
     *Find the direct neighbors of each vertex
     */
    protected void needNeighbors(){
	
	if (neighbors!= null && neighbors.length == vertices.size())
	    return;

	int vidx;

	neighbors= new Vector[vertices.size()];
	for (int i = 0; i < vertices.size(); i++)
	    neighbors[i]= new Vector();

	for (int s = 0; s < modelSegments.size(); s++){
	    ModelSegment tmpModelSegment = (ModelSegment) modelSegments.elementAt(s);
	    for (int f = 0; f < tmpModelSegment.elements.size(); f++) {
		ModelElement tmpf = (ModelElement) (tmpModelSegment.elements.elementAt(f));					
		if (tmpf.indices.size() > 0){ 
		    for (int fp = 0; fp < 3; fp++){
			vidx = tmpf.indices.elementAt(fp).intValue();
			Vector<Integer> me= neighbors[vidx-1];
			Integer n1 = tmpf.indices.elementAt((fp+1)%3);
			Integer n2 = tmpf.indices.elementAt((fp+2)%3);
			if (me.indexOf(n1)==-1)
			    me.add(n1);
			if (me.indexOf(n2)==-1)
			    me.add(n2);
		    }
		}
	    }
	}
    }
    
    
    /**
     *Find the faces touching each vertex
     */
    protected void needAdjacentFaces(){

	if (adjacentFaces!=null && adjacentFaces.length==vertices.size())
	    return;
	
        int vidx;
	
	adjacentFaces= new Vector[vertices.size()];
	for (int i = 0; i < vertices.size(); i++){
	    adjacentFaces[i]= new Vector();
	}
	
	for (int s = 0; s < modelSegments.size(); s++){
	    ModelSegment tmpModelSegment = (ModelSegment) modelSegments.elementAt(s);
	    for (int f = 0; f < tmpModelSegment.elements.size(); f++) {
		ModelElement tmpf = (ModelElement) (tmpModelSegment.elements.elementAt(f));					
		if (tmpf.indices.size() > 0){ 
		    for (int fp = 0; fp < 3; fp++){
			vidx = tmpf.indices.elementAt(fp).intValue();
			adjacentFaces[vidx-1].add(new SegmentAndFace(s, f));
		    }
		}
	    }
	}
    }
    
    

    protected class SegmentAndFace{

	int segment;
	int face;

	SegmentAndFace(int sf){
	    segment= sf;
	    face= sf;
	}

	SegmentAndFace(int s, int f){
	    segment= s;
	    face= f;
	}
	
	public String toString(){
	    return "segment: " + segment + " face: " + face;
	}

	public boolean equals(Object obj){
	    return segment==((SegmentAndFace)obj).segment &&
		face==((SegmentAndFace)obj).face;
	} 
    }

    /**
     * Find the face across each edge from each other face (-1 on boundary)
     * If topology is bad, not necessarily what one would expect...
     */
    protected void needAcrossEdge(){
	if (acrossEdge!=null && acrossEdge.size()==modelSegments.size())
	    return;
	
	needAdjacentFaces();
	
	acrossEdge= new Vector();
	for(int s= 0; s<modelSegments.size(); s++){
	    ModelSegment tmpModelSegment= modelSegments.elementAt(s);
	    acrossEdge.add(new SegmentAndFace[tmpModelSegment.elements.size()][3]);
	    
	    for(int f = 0; f < tmpModelSegment.elements.size(); f++) {
		((SegmentAndFace[][])acrossEdge.elementAt(s))[f][0]= new SegmentAndFace(-1);
		((SegmentAndFace[][])acrossEdge.elementAt(s))[f][1]= new SegmentAndFace(-1);
		((SegmentAndFace[][])acrossEdge.elementAt(s))[f][2]= new SegmentAndFace(-1);
	    }
	}
	
	
	for (int s = 0; s < modelSegments.size(); s++){
	    ModelSegment tmpModelSegment = (ModelSegment) modelSegments.elementAt(s);
	    for (int f = 0; f < tmpModelSegment.elements.size(); f++) {
		ModelElement tmpf = (ModelElement) (tmpModelSegment.elements.elementAt(f));
		for (int fp = 0; fp < 3; fp++){
		    if (((SegmentAndFace[][])acrossEdge.elementAt(s))[f][fp].face != -1)
			continue;
		    int v1 = tmpf.indices.elementAt((fp+1)%3).intValue();
		    int v2 = tmpf.indices.elementAt((fp+2)%3).intValue();
		    Vector<SegmentAndFace> a1 = adjacentFaces[v1-1];
		    Vector<SegmentAndFace> a2 = adjacentFaces[v2-1];
		    for (int k1 = 0; k1 < a1.size(); k1++) {
			SegmentAndFace other= a1.elementAt(k1);
			if (other.segment == s && other.face == f)
			    continue;
			if (a2.indexOf(other) == -1)
			continue;
			ModelSegment segment= modelSegments.elementAt(other.segment);
			ModelElement face= segment.elements.elementAt(other.face);
			int ind= (face.indices.indexOf(v1)+1)%3;
			if(face.indices.elementAt((ind+1)%3).intValue() != v2)
			    continue;
			((SegmentAndFace[][])acrossEdge.elementAt(s))[f][fp].segment= other.segment;
			((SegmentAndFace[][])acrossEdge.elementAt(s))[f][fp].face= other.face;
			((SegmentAndFace[][])acrossEdge.elementAt(other.segment))[other.face][ind].segment= s;
			((SegmentAndFace[][])acrossEdge.elementAt(other.segment))[other.face][ind].face= f;
			break;
		    }
		}
	    }
	}
    }  


    /**
     * Compute per-vertex point areas
     */
    protected void needPointareas(){
	if (pointareas!= null && pointareas.length == vertices.size())
	    return;

	pointareas= new float[vertices.size()];
      	cornerareas= new Vector(modelSegments.size());
	for(int s= 0; s<modelSegments.size(); s++){
	    ModelSegment tmpModelSegment = (ModelSegment) modelSegments.elementAt(s);
	    cornerareas.add(new float[tmpModelSegment.elements.size()][3]);
	}
	
	for (int s = 0; s < modelSegments.size(); s++){
	    ModelSegment tmpModelSegment = (ModelSegment) modelSegments.elementAt(s);
	    for (int f = 0; f < tmpModelSegment.elements.size(); f++) {
		ModelElement tmpf = (ModelElement) (tmpModelSegment.elements.elementAt(f));
		int v0 = (tmpf.indices.elementAt(0)).intValue();
		int v1 = (tmpf.indices.elementAt(1)).intValue();
		int v2 = (tmpf.indices.elementAt(2)).intValue();
		
		
		PVector e[] = { PVector.sub(vertices.elementAt(v2-1), vertices.elementAt(v1-1)),
				PVector.sub(vertices.elementAt(v0-1), vertices.elementAt(v2-1)),
				PVector.sub(vertices.elementAt(v1-1), vertices.elementAt(v0-1)) };
				
		
		// Compute corner weights
		float area = 0.5f * (e[1].cross(e[0])).mag();
		float l2[] = { parent.sq((e[0]).mag()), 
			       parent.sq((e[1]).mag()), 
			       parent.sq((e[2]).mag()) };
		float ew[] = { l2[0] * (l2[1] + l2[2] - l2[0]),
			       l2[1] * (l2[2] + l2[0] - l2[1]),
			       l2[2] * (l2[0] + l2[1] - l2[2]) };
		if (ew[0] <= 0.0f) {
		    cornerareas.elementAt(s)[f][1]= -0.25f * l2[2] * area /
			e[0].dot(e[2]);
		    cornerareas.elementAt(s)[f][2]= -0.25f * l2[1] * area /
			e[0].dot(e[1]);
		    cornerareas.elementAt(s)[f][0]= area - 
			cornerareas.elementAt(s)[f][1] - 
			cornerareas.elementAt(s)[f][2];
		} else if (ew[1] <= 0.0f) {
		    cornerareas.elementAt(s)[f][2]= -0.25f * l2[0] * area /
			(e[1].dot(e[0]));
		    cornerareas.elementAt(s)[f][0]= -0.25f * l2[2] * area /
			(e[1].dot(e[2]));
		    cornerareas.elementAt(s)[f][1]= area - 
			cornerareas.elementAt(s)[f][2] - 
			cornerareas.elementAt(s)[f][0];
		} else if (ew[2] <= 0.0f) {
		    cornerareas.elementAt(s)[f][0]= -0.25f * l2[1] * area /
			(e[2].dot(e[1]));
		    cornerareas.elementAt(s)[f][1]= -0.25f * l2[0] * area /
			(e[2].dot(e[0]));
		    cornerareas.elementAt(s)[f][2]= area - 
			cornerareas.elementAt(s)[f][0] - 
			cornerareas.elementAt(s)[f][1];
		} else {
		    float ewscale = 0.5f * area / (ew[0] + ew[1] + ew[2]);
		    for (int j = 0; j < 3; j++)
			cornerareas.elementAt(s)[f][j] = ewscale * (ew[(j+1)%3] + ew[(j+2)%3]);
		}
		pointareas[v0-1] += cornerareas.elementAt(s)[f][0];
		pointareas[v1-1] += cornerareas.elementAt(s)[f][1];
		pointareas[v2-1] += cornerareas.elementAt(s)[f][2];
	    }
	}
    }

    /**
     * Rotate a coordinate system to be perpendicular to the given normal
     */
    protected void rotCoordSys(PVector old_u, PVector old_v,
			    PVector newNorm,
			    PVector new_u, PVector new_v){
	
	new_u.x = old_u.x;
	new_u.y = old_u.y;
	new_u.z = old_u.z;
	new_v.x = old_v.x;
	new_v.y = old_v.y;
	new_v.z = old_v.z;
	PVector oldNorm = old_v.cross(old_u);
	float ndot = oldNorm.dot(newNorm);
	if (ndot <= -1.0f) {
	    new_u.mult(-1);
	    new_v.mult(-1);
	    return;
	}
	PVector perpOld = PVector.sub(newNorm, PVector.mult(oldNorm, ndot));
	//PVector dperp = PVector.div(new PVector(1.0f, 1.0f, 1.0f), PVector.mult(PVector.add(oldNorm, newNorm), (1 + ndot)));
	PVector dperp = PVector.mult(PVector.add(oldNorm, newNorm), 1.0f/(1 + ndot));
	new_u.sub(PVector.mult(dperp, (new_u.dot(perpOld))));
	new_v.sub(PVector.mult(dperp, (new_v.dot(perpOld))));
    }


    /* Reproject a curvature tensor from the basis spanned by old_u and old_v (which are assumed to be unit-length and perpendicular) to the new_u, new_v basis.
     */
    protected void projCurv(PVector old_u, PVector old_v,
		  float old_ku, float old_kuv, float old_kv,
		  PVector new_u, PVector new_v,
		  PVector new_kuv){

	PVector r_new_u= new PVector(), r_new_v= new PVector();
	rotCoordSys(new_u, new_v, old_v.cross(old_u), r_new_u, r_new_v);
	float u1 = r_new_u.dot(old_u);
	float v1 = r_new_u.dot(old_v);
	float u2 = r_new_v.dot(old_u);
	float v2 = r_new_v.dot(old_v);
	
	new_kuv.x= old_ku * u1*u1 + old_kuv * (2.0f  * u1*v1) + old_kv * v1*v1;
	new_kuv.y= old_ku * u1*u2 + old_kuv * (u1*v2 + u2*v1) + old_kv * v1*v2;
	new_kuv.z= old_ku * u2*u2 + old_kuv * (2.0f  * u2*v2) + old_kv * v2*v2;
    }



    /**
     * Like projCurv, but for dcurv
     */
    protected void projDcurv(PVector old_u, PVector old_v,
		   float old_dcurv[],
		   PVector new_u, PVector new_v,
		   float new_dcurv[]){
	
	PVector r_new_u= new PVector(), r_new_v= new PVector();
	rotCoordSys(new_u, new_v, old_v.cross(old_u), r_new_u, r_new_v);

	float u1 = r_new_u.dot(old_u);
	float v1 = r_new_u.dot(old_v);
	float u2 = r_new_v.dot(old_u);
	float v2 = r_new_v.dot(old_v);

	new_dcurv[0] = old_dcurv[0]*u1*u1*u1 +
		       old_dcurv[1]*3.0f*u1*u1*v1 +
		       old_dcurv[2]*3.0f*u1*v1*v1 +
		       old_dcurv[3]*v1*v1*v1;
	new_dcurv[1] = old_dcurv[0]*u1*u1*u2 +
		       old_dcurv[1]*(u1*u1*v2 + 2.0f*u2*u1*v1) +
		       old_dcurv[2]*(u2*v1*v1 + 2.0f*u1*v1*v2) +
		       old_dcurv[3]*v1*v1*v2;
	new_dcurv[2] = old_dcurv[0]*u1*u2*u2 +
		       old_dcurv[1]*(u2*u2*v1 + 2.0f*u1*u2*v2) +
		       old_dcurv[2]*(u1*v2*v2 + 2.0f*u2*v2*v1) +
		       old_dcurv[3]*v1*v2*v2;
	new_dcurv[3] = old_dcurv[0]*u2*u2*u2 +
		       old_dcurv[1]*3.0f*u2*u2*v2 +
		       old_dcurv[2]*3.0f*u2*v2*v2 +
		       old_dcurv[3]*v2*v2*v2;
    }



    /**
     * Given a curvature tensor, find principal directions and curvatures
     * Makes sure that pdir1 and pdir2 are perpendicular to normal
     */
    protected void diagonalizeCurv(PVector old_u, PVector old_v,
			 float ku, float kuv, float kv,
			 PVector new_norm,
			 PVector pdir1, PVector pdir2, float k[]){
	
	PVector pdir1new, pdir2new;

	PVector r_old_u= new PVector(), r_old_v= new PVector();
	rotCoordSys(old_u, old_v, new_norm, r_old_u, r_old_v);
	float c = 1, s = 0, tt = 0;
	if (kuv != 0.0f) {
	    // Jacobi rotation to diagonalize
	    float h = 0.5f * (kv - ku) / kuv;
	    tt = (h < 0.0f) ?
		1.0f / (h - (float)Math.sqrt(1.0f + h*h)) :
		1.0f / (h + (float)Math.sqrt(1.0f + h*h));
	    c = 1.0f / (float)Math.sqrt(1.0f + tt*tt);
	    s = tt * c;
	}

	k[0] = ku - tt * kuv;
	k[1] = kv + tt * kuv;

	if (Math.abs(k[0]) >= Math.abs(k[1])) {
	    pdir1new = PVector.sub(PVector.mult(r_old_u, c),  PVector.mult(r_old_v, s));
	} 
	else {
	    float k0= k[0];
	    k[0]= k[1];
	    k[1]= k0;
	    pdir1new = PVector.add(PVector.mult(r_old_u, s),  PVector.mult(r_old_v, c));
	}
	pdir2new = pdir1new.cross(new_norm);
	pdir1.x= pdir1new.x;
	pdir1.y= pdir1new.y;
	pdir1.z= pdir1new.z;
	pdir2.x= pdir2new.x;
	pdir2.y= pdir2new.y;
	pdir2.z= pdir2new.z;
    }


    /**
     * Perform LDL^T decomposition of a symmetric positive definite matrix.
     * Like Cholesky, but no square roots.  Overwrites lower triangle of matrix.
     */
    protected boolean ldltdc3(float A[][], float rdiag[]){
	float v[]= new float[2];
	for (int i = 0; i < 3; i++) {
	    for (int k = 0; k < i; k++)
		v[k] = A[i][k] * rdiag[k];
	    for (int j = i; j < 3; j++) {
		float sum = A[i][j];
		for (int k = 0; k < i; k++)
		    sum -= v[k] * A[j][k];
		if (i == j) {
		    if (sum <= 0)
			return false;
		    rdiag[i] = 1 / sum;
		} else {
		    A[j][i] = sum;
		}
	    }
	}
	
	return true;
    }


    /**
     * Perform LDL^T decomposition of a symmetric positive definite matrix.
     * Like Cholesky, but no square roots.  Overwrites lower triangle of matrix.
     */
    protected boolean ldltdc4(float A[][], float rdiag[]){
	float v[]= new float[3];
	for (int i = 0; i < 4; i++) {
	    for (int k = 0; k < i; k++)
		v[k] = A[i][k] * rdiag[k];
	    for (int j = i; j < 4; j++) {
		float sum = A[i][j];
		for (int k = 0; k < i; k++)
		    sum -= v[k] * A[j][k];
		if (i == j) {
		    if (sum <= 0)
			return false;
		    rdiag[i] = 1 / sum;
		} else {
		    A[j][i] = sum;
		}
	    }
	}
	
	return true;
    }


    /** 
     * Solve Ax=B after ldltdc
     */
    protected void ldltsl3(float A[][], float rdiag[], float B[], float x[]){
	
	int i;
	for (i= 0; i < 3; i++) {
	    float sum= B[i];
	    for (int k= 0; k < i; k++)
		sum-= A[i][k] * x[k];
	    x[i]= sum * rdiag[i];
	}
	for (i= 2; i >= 0; i--) {
	    float sum= 0;
	    for (int k= i + 1; k < 3; k++)
		sum+= A[k][i] * x[k];
	    x[i]-= sum * rdiag[i];
	}
    }
    



    /** 
     * Solve Ax=B after ldltdc
     */
    protected void ldltsl4(float A[][], float rdiag[], float B[], float x[]){
	
	int i;
	for (i = 0; i < 4; i++) {
	    float sum= B[i];
	    for (int k= 0; k < i; k++)
		sum -= A[i][k] * x[k];
	    x[i] = sum * rdiag[i];
	}
	for (i = 3; i >= 0; i--) {
	    float sum= 0;
	    for (int k= i + 1; k < 4; k++)
		sum+= A[k][i] * x[k];
	    x[i]-= sum * rdiag[i];
	}
    }
    

    /**
     * Compute principal curvatures and directions.
     */
    protected void needCurvatures(){
	if (curv1!= null && curv1.length == vertices.size())
	    return;
	
	if(!isTriangleMesh) return;

	needPointareas();

	int nv= vertices.size();

	// Resize the arrays we'll be using
	curv1= new float[nv];
	curv2= new float[nv];
	pdir1= new PVector[nv];
	pdir2= new PVector[nv];
	for(int i=0; i<nv; i++){
	    pdir1[i]= new PVector();
	    pdir2[i]= new PVector();
	}
	float curv12[]= new float[nv];

	PVector normals[]= new PVector[nv];

	// Set up an initial coordinate system per vertex
	for (int s = 0; s < modelSegments.size(); s++){
	    ModelSegment tmpModelSegment = modelSegments.elementAt(s);
	    for (int f = 0; f < tmpModelSegment.elements.size(); f++) {
		ModelElement tmpf = tmpModelSegment.elements.elementAt(f);
		int v0 = (tmpf.indices.elementAt(0)).intValue();
		int v1 = (tmpf.indices.elementAt(1)).intValue();
		int v2 = (tmpf.indices.elementAt(2)).intValue();
		
		
		pdir1[v0-1]= PVector.sub(vertices.elementAt(v1-1),
					 vertices.elementAt(v0-1));
		pdir1[v1-1]= PVector.sub(vertices.elementAt(v2-1), 
					 vertices.elementAt(v1-1));
		pdir1[v2-1]= PVector.sub(vertices.elementAt(v0-1), 
					 vertices.elementAt(v2-1));

		for(int i= 0; i<3; i++){
		    int v = (tmpf.indices.elementAt(i)).intValue();
		    int n = (tmpf.nindices.elementAt(i)).intValue();
		    normals[v-1]= normv.elementAt(n-1);
		} 
	    }
	}
	
	for(int i= 0; i<nv; i++){
	    pdir1[i]= normals[i].cross(pdir1[i]);
	    pdir1[i].normalize();
	    pdir2[i]= pdir1[i].cross(normals[i]);
	}
	
	// Compute curvature per-face
	for(int s = 0; s < modelSegments.size(); s++){
	    ModelSegment tmpModelSegment = modelSegments.elementAt(s);
	    for (int f = 0; f < tmpModelSegment.elements.size(); f++) {
		ModelElement tmpf = tmpModelSegment.elements.elementAt(f);
	    
		int v0 = (tmpf.indices.elementAt(0)).intValue();
		int v1 = (tmpf.indices.elementAt(1)).intValue();
		int v2 = (tmpf.indices.elementAt(2)).intValue();
	    
		// Edges
		PVector e[] = { PVector.sub(vertices.elementAt(v2-1),
					    vertices.elementAt(v1-1)),
				PVector.sub(vertices.elementAt(v0-1),
					    vertices.elementAt(v2-1)),
				PVector.sub(vertices.elementAt(v1-1),
					    vertices.elementAt(v0-1)) };
	    
		// N-T-B coordinate system per face
		PVector t = e[0].get();
		t.normalize();
		PVector n = e[1].cross(e[0]);
		PVector b = t.cross(n);
		b.normalize();
	    
		// Estimate curvature based on variation of normals
		// along edges
		float m[] = { 0, 0, 0 };
		float w[][] = { {0,0,0}, {0,0,0}, {0,0,0} };
		for (int j = 0; j < 3; j++) {
		    float u = e[j].dot(t);
		    float v = e[j].dot(b);
		    w[0][0] += u*u;
		    w[0][1] += u*v;
		    //w[1][1] += v*v + u*u; 
		    //w[1][2] += u*v; 
		    w[2][2] += v*v;
		    int prevj= j>0 ? j-1 : j+2;
		    PVector dn = PVector.sub(normv.elementAt(tmpf.nindices.elementAt(prevj).intValue()-1),
					     normv.elementAt(tmpf.nindices.elementAt((j+1)%3).intValue()-1));
		    float dnu = dn.dot(t);
		    float dnv = dn.dot(b);
		    m[0] += dnu*u;
		    m[1] += dnu*v + dnv*u;
		    m[2] += dnv*v;
	
		}
		w[1][1] = w[0][0] + w[2][2];
		w[1][2] = w[0][1];
	    
		// Least squares solution
		float diag[]= new float[3];
		if (!ldltdc3(w, diag)) {
		    //System.out.println("ldltdc failed!");
		    continue;
		}
		ldltsl3(w, diag, m, m);
		
		// Push it back out to the vertices
		for (int j = 0; j < 3; j++) {
		    int vj = tmpf.indices.elementAt(j).intValue();
		    PVector c= new PVector();
		    projCurv(t, b, m[0], m[1], m[2],
			     pdir1[vj-1], pdir2[vj-1], c);

		    float wt = cornerareas.elementAt(s)[f][j] / pointareas[vj-1];
		    curv1[vj-1]  += wt * c.x;
		    curv12[vj-1] += wt * c.y;
		    curv2[vj-1]  += wt * c.z;
		}
	    }
	}
	float curvi[]= new float[2];
	// Compute principal directions and curvatures at each vertex

	for(int i= 0; i< nv; i++){
	    diagonalizeCurv(pdir1[i], pdir2[i],
			    curv1[i], curv12[i], curv2[i],
			    normals[i], pdir1[i], pdir2[i],
			    curvi);
	    curv1[i]= curvi[0];
	    curv2[i]= curvi[1];
	}
    }


    // Compute derivatives of curvature.
    protected void needDCurv(){
	
	if (dcurv!= null && dcurv.length == vertices.size())
	    return;
	needCurvatures();
	
	dcurv= new float[vertices.size()][4];
	
	// Compute dcurv per-face
	for(int s = 0; s < modelSegments.size(); s++){
	    ModelSegment tmpModelSegment = (ModelSegment) modelSegments.elementAt(s);
	    for (int f = 0; f < tmpModelSegment.elements.size(); f++) {
		ModelElement tmpf = (ModelElement) (tmpModelSegment.elements.elementAt(f));
		
		int vidx0 = (tmpf.indices.elementAt(0)).intValue();
		int vidx1 = (tmpf.indices.elementAt(1)).intValue();
		int vidx2 = (tmpf.indices.elementAt(2)).intValue();
		
		// Edges
		PVector e[] = { PVector.sub(vertices.elementAt(vidx2-1),
					    vertices.elementAt(vidx1-1)),
				PVector.sub(vertices.elementAt(vidx0-1),
					    vertices.elementAt(vidx2-1)),
				PVector.sub(vertices.elementAt(vidx1-1),
					    vertices.elementAt(vidx0-1)) };

		// N-T-B coordinate system per face
		PVector t = e[0].get();
		t.normalize();
		PVector n = e[1].cross(e[0]);
		PVector b = t.cross(n);
		b.normalize();
		
		// Project curvature tensor from each vertex into this
		// face's coordinate system
		PVector fcurv[]= new PVector[3]; 
		for (int j = 0; j < 3; j++) {
		    fcurv[j]= new PVector();
		    int vj = tmpf.indices.elementAt(j).intValue();
		    projCurv(pdir1[vj-1], pdir2[vj-1], curv1[vj-1], 0, curv2[vj-1], t, b, fcurv[j]);

		}

		// Estimate dcurv based on variation of curvature along edges
		float m[] = { 0, 0, 0, 0 };
		float w[][] = { {0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} };
		for (int j = 0; j < 3; j++) {
			// Variation of curvature along each edge
		    int prevfcurv= j>0 ? j-1 : j+2;
		    PVector dfcurv = PVector.sub(fcurv[prevfcurv], fcurv[(j+1)%3]);
		    float u = e[j].dot(t);
		    float v = e[j].dot(b);
		    float u2 = u*u, v2 = v*v, uv = u*v;
		    w[0][0] += u2;
		    w[0][1] += uv;
		    //w[1][1] += 2.0f*u2 + v2;
		    //w[1][2] += 2.0f*uv;
		    //w[2][2] += u2 + 2.0f*v2;
		    //w[2][3] += uv;
		    w[3][3] += v2;
		    m[0] += u*dfcurv.x;
		    m[1] += v*dfcurv.x + 2.0f*u*dfcurv.y;
		    m[2] += 2.0f*v*dfcurv.y + u*dfcurv.z;
		    m[3] += v*dfcurv.z;
		}
		w[1][1] = 2.0f * w[0][0] + w[3][3];
		w[1][2] = 2.0f * w[0][1];
		w[2][2] = w[0][0] + 2.0f * w[3][3];
		w[2][3] = w[0][1];

		// Least squares solution
		float d[]= new float[4];
		if (!ldltdc4(w, d)) {
		    //System.out.println("ldltdc failed!");
		    continue;
		}
		ldltsl4(w, d, m, m);
		float[] faceDcurv= new float[4];
		faceDcurv[0]= m[0];
		faceDcurv[1]= m[1];
		faceDcurv[2]= m[2];
		faceDcurv[3]= m[3];

		// Push it back out to each vertex
		for (int j = 0; j < 3; j++) {
		    int vj = tmpf.indices.elementAt(j).intValue();
			float thisVertDcurv[]= new float[4];
			projDcurv(t, b, faceDcurv,
				   pdir1[vj-1], pdir2[vj-1], thisVertDcurv);
			float wt = cornerareas.elementAt(s)[f][j] / pointareas[vj-1];
			dcurv[vj-1][0]+= wt * thisVertDcurv[0];
			dcurv[vj-1][1]+= wt * thisVertDcurv[1];
			dcurv[vj-1][2]+= wt * thisVertDcurv[2];
			dcurv[vj-1][3]+= wt * thisVertDcurv[3];
		}
	    }
	}
    }



    /**
     * Convert faces to tstrips
     */
    protected void needTStrips(){

	if (tstrips!= null && !tstrips.isEmpty())
		return;

	needAcrossEdge();

	Vector<SegmentAndFace> todo= new Vector<SegmentAndFace>();
	int numFaces= 0;
	Vector<SegmentAndFace> faces= new Vector<SegmentAndFace>();
	Vector<byte[]>faceAvail= new Vector<byte[]>(modelSegments.size());

	for(int s= 0; s<modelSegments.size(); s++){
	    ModelSegment tmpModelSegment = (ModelSegment) modelSegments.elementAt(s);
	    numFaces+= tmpModelSegment.elements.size();
	    for (int f = 0; f < tmpModelSegment.elements.size(); f++) {
		faces.add(new SegmentAndFace(s, f));
	    }
	    faceAvail.add(new byte[tmpModelSegment.elements.size()]);
	}

	for(int s = 0; s < modelSegments.size(); s++){
	    ModelSegment tmpModelSegment = (ModelSegment) modelSegments.elementAt(s);
	    for (int f = 0; f < tmpModelSegment.elements.size(); f++) {
		
		faceAvail.elementAt(s)[f] = 
		    (byte)(((acrossEdge.elementAt(s)[f][0].face != -1) ? 1 : 0) +
			   ((acrossEdge.elementAt(s)[f][1].face != -1) ? 1 : 0) +
			   ((acrossEdge.elementAt(s)[f][2].face != -1) ? 1 : 0));
		if (faceAvail.elementAt(s)[f] == 1)
			todo.add(new SegmentAndFace(s, f));
	    }
	}
	
	tstrips= new Vector<Integer>(numFaces*2);

	int nstrips = 0;
	int f = 0;
	while (f < numFaces) {
	    SegmentAndFace next;
	    if (todo.isEmpty()) {
		next= faces.elementAt(f++);
	    } else {
		next = todo.lastElement();
		todo.remove(todo.size()-1);
	    }
	    if (faceAvail.elementAt(next.segment)[next.face] < 0)
		continue;
	    tstripBuild(next, faceAvail, todo);
	    nstrips++;
	}
	
	convertStrips(TStripRep.TSTRIP_LENGTH);
    }
    

    /**
     * Build a triangle strip starting with the given face
     */
     protected void tstripBuild(SegmentAndFace sf, Vector<byte[]> faceAvail,
	Vector<SegmentAndFace> todo){
	 int s= sf.segment, f= sf.face;
	 ModelElement face= modelSegments.elementAt(s).elements.elementAt(f);
	 if (faceAvail.elementAt(s)[f] == 0){
	     tstrips.add(face.indices.elementAt(0));
	     tstrips.add(face.indices.elementAt(1));
	     tstrips.add(face.indices.elementAt(2));
	     tstrips.add(Integer.valueOf(-1));
	     faceAvail.elementAt(s)[f] = -1;
	     return;
	 }

	 int score[]= new int[3];
	 for (int i = 0; i < 3; i++) {
	     score[i] = 0;
	     SegmentAndFace ae = acrossEdge.elementAt(s)[f][i];
	     if (ae.face == -1 || faceAvail.elementAt(ae.segment)[ae.face] < 0)
		 continue;
	     score[i]++;
	     int nextEdge = modelSegments.elementAt(ae.segment).elements.elementAt(ae.face).indices.indexOf(face.indices.elementAt((i+1)%3));
	     SegmentAndFace nae = acrossEdge.elementAt(ae.segment)[ae.face][nextEdge];
	     if (nae.face == -1 || faceAvail.elementAt(nae.segment)[nae.face] < 0)
		 continue;
	     score[i]++;
	     if (faceAvail.elementAt(ae.segment)[ae.face] == 2)
		 score[i]++;
	 }
	 
	 int best_score = Math.max(Math.max(score[0], score[1]), score[2]);
	 int best = (score[0] == best_score) ? 0 :
	     (score[1] == best_score) ? 1 : 2;
	 
	 Integer vlast2 = face.indices.elementAt( best     );
	 Integer vlast1 = face.indices.elementAt((best+1)%3);
	 Integer vnext  = face.indices.elementAt((best+2)%3);
	 int dir = 1;
	 tstrips.add(vlast2);
	 tstrips.add(vlast1);
	 
	 while (true) {
	     tstrips.add(vnext);
	     faceAvail.elementAt(s)[f] = -1;
		for (int j = 0; j < 3; j++) {
		    SegmentAndFace ae = acrossEdge.elementAt(s)[f][j];
		    if (ae.face == -1)
			continue;
		    if (faceAvail.elementAt(ae.segment)[ae.face] > 0)
			faceAvail.elementAt(ae.segment)[ae.face]--;
		    if (faceAvail.elementAt(ae.segment)[ae.face] == 1)
				todo.add(ae);
		}
	    
		s = acrossEdge.elementAt(s)[f][modelSegments.elementAt(s).elements.elementAt(f).indices.indexOf(vlast2)].segment;
		f = acrossEdge.elementAt(sf.segment)[f][modelSegments.elementAt(sf.segment).elements.elementAt(f).indices.indexOf(vlast2)].face;
		if (f == -1 || faceAvail.elementAt(s)[f] < 0)
		    break;
		vlast2 = vlast1;
		vlast1 = vnext;
		vnext = modelSegments.elementAt(s).elements.elementAt(f).indices.elementAt((modelSegments.elementAt(s).elements.elementAt(f).indices.indexOf(vlast2)+3+dir)%3);
		dir = -dir;
	 }
	 
	 tstrips.add(Integer.valueOf(-1));
     }
    

    /**
     * Convert between "length preceding strip" and "-1 following strip" representations
     */
    protected void convertStrips(TStripRep rep){
    
	if (tstrips.isEmpty())
		return;
	if (rep == TStripRep.TSTRIP_TERM && tstrips.lastElement().intValue() == -1)
		return;
	if (rep == TStripRep.TSTRIP_LENGTH && tstrips.lastElement().intValue() != -1) {
		return;
	}

	if (rep == TStripRep.TSTRIP_TERM) {
	    int len = tstrips.get(0).intValue();
	    for (int i = 1; i < tstrips.size(); i++) {
		if (len>0) {
		    tstrips.set(i-1, tstrips.get(i));
		    len--;
		} else {
		    tstrips.set(i-1, Integer.valueOf(-1));
		    len = tstrips.get(i);
		}
	    }
	    tstrips.set(tstrips.size()-1, Integer.valueOf(-1));
	} 
	else {
	    int len = 0;
	    for (int i = tstrips.size() - 2; i >= 0; i--){
		if (tstrips.get(i).intValue() == -1){
		    tstrips.set(i+1, Integer.valueOf(len));
		    len = 0;
		} 
		else {
		    tstrips.set(i+1, tstrips.get(i));
		    len++;
		}
	    }
	    tstrips.set(0, Integer.valueOf(len));
	    //collect_tris_in_strips(tstrips);
	}
    }
}
